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I.  Introduction 

The  development  of  biomedical  optical  imaging  for  the  detection  of  breast  cancer 
has  been  the  subject  of  an  increasing  number  of  research  programs  in  science, 
engineering,  and  medicine.  The  ability  to  optically  detect  breast  cancer  in  a  group  of 
women  who  are  not  good  candidates  for  conventional  mammography  depends  upon  the 
optical  contrast  available  between  normal  and  diseased  tissues.  In  a  study  published  in 
1996,  we  measured  the  optical  properties  of  115  human  breast  tissue  specimens  in  order 
to  ascertain  the  level  of  contrast  which  may  be  available  for  the  proposed  general 
screening  tool  of  “optical  mammography.”  Our  results,  although  they  underpredict  the 
contrast  owing  to  absorption,  suggest  that  there  is  insufficient  contrast  for  the  detection  of 
disease.  In  addition,  that  study  suggests  the  need  for  an  optical  contrast  agent. 

The  detection  of  superficially  located  diseased  tissues  via  the  use  of  fluorescent 
dyes  and  photodynamic  drugs  as  contrast  agents  has  been  proposed  by  several 
investigators  over  the  past  several  years.  In  our  program,  we  have  investigated  the 
development  of  imaging  technologies  which  employ  fluorescent  drugs  or  agents  for 
detection  of  small  lesions  beneath  the  skin,  or  tissue  volumes  that  are  located  several 
centimeters  deep,  However,  to  accomplish  such  a  technology,  a  long-standing  problem 
must  be  overcome.  Specifically,  the  low  uptake  or  “leakage”  into  diseased  tissues, 
providing  insufficient  contrast  for  detection  of  disease,  While  targeted  delivery  may 
improve  uptake  ratios  and  contrast  imparted  by  these  agents,  drug  and  contrast  agent 
targeting  approaches  have  been  rather  elusive  even  in  conventional  imaging  modalities 
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such  as  magnetic  resonance  imaging  and  x-ray  computer  tomography.  Investigators  have 
sought  to  improve  optical  contrast  independent  of  uptake  ratio  by  employing  dyes  which 
fluoresce  differently  in  diseased  and  normal  tissues.  The  use  of  agents  whose  emission 
characteristics  vary  with  tissue  pH  and  02  may  not  only  enhance  detection  of  diseased 
tissues  by  the  nature  of  differing  fluorescent  properties,  but  may  also  contain  diagnostic 
information  regarding  the  disease. 

The  difficulty,  however,  lies  in  using  the  multiply  scattered  light  detected  at  the 
tissue-air  interface  to  reconstruct  and  image  which  differentiates  internal  tissues  on  the 
basis  of  uptake  and  the  fluorescent  properties  of  quantum  yield  and  lifetime.  In  ongoing 
work  in  this  and  other  laboratories,  investigators  have  demonstrated  the  capacity  for 
fluorescent  lifetime  and  quantum  yield  imaging  form  frequency-domain  measurements  of 
modulation  phase-shift  and  amplitude  demodulation  conducted  at  the  surface  of 
simulated  tissue  phantoms.  Under  this  support,  the  aims  of  the  exploratory  project  were: 

(1)  To  develop  inverse  imaging  algorithms  providing  lifetime  and  yield  mapping  of 
tissues  as  thick  as  five  centimeters; 

(2)  To  measure  the  temporal,  re-emitted  excitation  and  fluorescence  from  tissue 
phantoms  using  frequency-domain  techniques;  and 

(3)  To  use  these  measurements  as  input  into  the  inverse  imaging  algorithm  for 
reconstruction  of  lifetime  and  yield  maps. 
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Specifically,  we  seek  to  develop  fluorescence  lifetime  imaging  techniques  for 
screening  in  a  population  of  women  who  are  at  high  risk  for  the  disease,  and  for  those 
women  who  require  prognostic  information  of  lymphatic  involvement  in  order  to  aid  in 
the  management  of  the  disease.  This  DOD  support  focuses  upon  the  development  of 
inverse  solution  procedures  and  coupling  of  measurements  to  image  algorithms  for  the 
detection  of  disease.  These  specific  aims  are  the  foundations  of  the  technology. 

n.  Computational  Methods 

As  described  in  the  publications  in  the  Appendix,  our  work  has  been  towards  the 
development  of  forward  simulators  of  the  propagation  of  excitation  and  emission  light  in 
tissues.  Under  this  support,  we  have  generated  two  forward  solvers: 

(1)  a  finite  difference,  MUDPACK  routine  operating  in  two  and  three  dimensions,  and 

(2)  a  finite  element  approach  which  enables  arbitrary  tissue  geometry. 

We  have  also  developed  Newton-Raphson  based  inversion  strategies  which 
enable  the  reconstruction  of  fluorescent  lifetime  and  the  product  of  quantum  efficiency 
and  absorption  cross  section  (Paithankar,  et  al.,  1997)  as  well  as  the  reconstruction  of 
fluorescent  lifetime,  quantum  efficiency,  and  absorption  cross  section  (Troy  and  Sevick, 
1998).  The  publications  are  included  in  the  Appendix  for  further  information. 
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More  recently,  in  the  latter  two  months  of  the  DOD  project,  we  have  coupled  with 
University  of  Vermont’s  Computer  Science  Department  in  order  to  develop  Bayesian 
optical  tomography  techniques.  The  approach  provides  a  zonation  structure  to  the  inverse 
optical  imaging  problem,  enabling  rapid  convergence  and  stability  for  reconstruction  of 
absorption  and  fluorescent  properties. 

HI.  Experimental  Methods 

In  our  inversion  strategies,  we  have  employed  synthetic  data  to  demonstrate  the 
feasibility  of  optical  imaging  of  absorption  and  fluorescence  following  the  administration 
of  a  contrast  agent.  In  order  to  demonstrate  the  feasibility  of  the  technology,  we  need  to 
develop  fast  methods  for  data  acquisition  of  frequency-domain  photon  migration 
measurements,  demonstrate  the  ability  to  match  experimental  measurements  with  those 
predicted  by  our  forward  solvers,  and  create  real,  experimental  data  sets  for  incorporation 
into  our  solutions  to  the  inverse  optical  imaging  problem.  In  part  under  this  support,  we 
have  developed  a  multi-pixel  imaging  device  which  consists  of  a  gain  modulated  image 
intensified  camera  to  capture  measurements  of  phase  and  amplitude  modulation  of  re¬ 
emitted  light  from  the  surface  of  tissues  without  the  need  to  optically  couple  or  contact 
the  tissue. 

The  details  of  the  apparatus  are  found  in  a  publication  included  in  the  appendix. 


Briefly,  the  technique  depends  upon  illuminating  the  tissue  volume  of  interest  with  a 
modulated  light  source,  i.e.,  either  a  point  source  of  light  delivered  by  a  fiber  optic,  or  a 
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plane  wave  created  by  the  expansion  of  a  diode  beam.  The  light  which  propagates  from 
the  source  within  the  tissue  and  to  the  various  positions  on  the  surface  of  the  tissue  is 
collected  by  the  ICCD  system  which  is  focused  upon  the  surface  of  interest.  Incidentally, 
the  ICCD  system  enables  in  principle  512  x  512  pixel  measurements  on  the  surface  of  a 
tissue  whose  dimensional  are  set  by  a  focusing  lens.  Measurements  can  be  made  of 
excitation  light,  as  well  as  from  the  fluorescence  generated  by  a  contrast  agent.  We  report 
that  the  measurements  are  hundreds  of  times  faster  than  single-pixel  scanning 
demonstrated  by  other  groups  and  imaging  companies  and  point  out  that  our 
measurements  contain  less  noise  than  that  reported  for  single-pixel  systems.  Since  the 
propagation  of  measurement  noise  into  an  inverse  solution  severely  deteriorates  the 
solution  or  the  determined  optical  image,  we  expect  that  our  measurement  approach 

provides  the  best  means  for  coupling  real  experimental  data  to  the  solution  for  the  inverse 
imaging  problem. 

While  we  were  unable  to  couple  the  measurements  with  the  inversion  algorithms 
owing  to  the  short  duration  and  small  funding  level  of  this  exploratory  research  proposal, 
we  have  nonetheless  demonstrated  the  ability  to  quantify  fluorescence  lifetime  in  a  highly 
scattering,  uniform  medium.  While  the  manuscript  are  currently  under  preparation  and 
will  be  made  available  to  the  funding  agency  upon  their  completion,  the  following  briefly 
summarizes  our  approach  and  findings.  These  experimental  measurements  confirm  the 
ability  to  quantify  fluorescence  lifetime  in  scattering  media. 
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The  goals  of  the  non-imaging  experiments  were  to  determine  the  fluorescent 
optical  properties  of  a  dye  in  a  highly  scattering  media.  Specifically,  we  sought  to 
determine  the  lifetime  of  Indocyanine  Green  (currently  used  a  s  a  contrast  agent  in  our 
animal  studies  supported  under  separate  funding)  and  IR-140,  a  common  laser  dye, 
Indocyanine  Green  has  a  measured  lifetime  of  0.58  ns  as  measured  in  our  laboratory 
while  IR-140  has  a  lifetime  of  1 .19  ns.  Using  frequency-domain  measurements  of  photon 
migration,  we  sought  to  determine  the  lifetime  of  these  two  dyes  when  contained  in 
micromolar  concentrations  within  a  tissue-like  scattering  solution  of  0.5%  Intralipid. 

This  gives  us  confidence  that  the  inverse  optical  imaging  of  fluorescence  lifetime  can  be 
achieved,  We  found  that  in  order  to  effectively  perform  lifetime  spectroscopy  in 
scattering  media,  we  required  three  FDPM  measurements:  One  measurement  was 
required  at  the  excitation  wavelength  in  which  a  neutral  density  filter  was  simply  placed 
at  the  detector  in  order  to  determine  the  propagation  characteristics  of  the  excitation  light. 
Another  measurement  involved  measuring  the  propagation  of  an  incident  wave  that  was 
at  the  same  wavelength  as  the  generated  fluorescence.  We  term  this  measurement  the 
emission  FDPM  measurement.  Finally,  the  final  measurement  consisted  of  exciting  the 
sample  with  the  excitation  light  and  using  interference  filters  to  collect  the  generated 
fluorescence  which  propagated  to  the  surface.  From  these  measurements,  we  were  able 
to  determine  the  lifetime  of  the  dye  without  using  any  a  priori  information  regarding  the 
optical  properties  of  the  medium.  Our  results,  again  which  will  be  made  available  to  the 
DOD  in  the  coming  weeks,  verify  the  feasibility  of  fluorescence  lifetime  imaging  in 
tissues. 
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IV.  CONCLUSIONS 


The  exploratory  research  supported  by  the  DOD  provided  the  following  results: 

(i)  The  inverse  imaging  problem  for  fluorescence  lifetime  has  a  solution  as  demonstrated 
using  synthetic  data  sets. 

(ii)  Multi-pixel  imaging  devices  were  developed  in  part  under  this  support  to  provide 
experimental  data  sets  for  demonstration  of  fluorescence  lifetime  imaging. 

(iii)  Using  single-pixel  measurements,  we  demonstrated  the  ability  to  perform  lifetime 
spectroscopy  in  homogeneous  scattering  media. 

In  the  future,  we  require  the  coupling  of  the  multi-pixel  measurements  with  the  solutions 
to  the  optical  imaging  modality  to  provide  a  new  imaging  modality  to  provide  detection 
and  perhaps  diagnostics  based  upon  spectroscopic  parameters  of  fluorescent 
photophysical  properties. 
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The  ability  to  map  interior  optical  properties  of  a  highly  scattering  medium  from 
exterior  measurements  of  light  propagation  is  afforded  by  optical  tomography.  In  this 
communication,  we  describe  the  problem  of  optical  tomography,  the  techniques  of 
photon  migration  measurements  necessary  to  accomplish  it,  and  the  development  of 
multipixel  measurements  for  rapid  collection  of  optical  signals.  These  multipixel 
measurements  are  shown  to  provide  detection  of  contrast  owing  to  the  optical  properties 
of  absorption  and  fluorescence  associated  with  dye-laden  heterogeneities  embedded 
in  a  tissue-like  scattering  medium.  From  these  rapid  measurements,  successful 
reconstruction  of  an  interior  optical  property  map  may  now  be  possible  with  clinically 
realistic  data  acquisition  times.  Applications  for  the  technology  arise  for  biomedical 
optical  imaging  for  the  in  vivo  detection  of  disease  and  the  diagnosis  of  tissue  (bio-) 
chemistry. 


Introduction 

Electrical  impedance,  ultrasound,  microwave,  and  opti¬ 
cal  tomographies  involve  launching  one  or  more  excita¬ 
tion  signals  at  the  periphery  of  a  body  whose  interior 
contents  are  under  investigation  and  measuring  the 
emitted  signals  at  a  number  of  detectors  or  receivers  also 
located  on  the  periphery.  Owing  to  the  presence  of 
interior  heterogeneities,  the  spatial  propagation  of  cur¬ 
rent,  sound,  microwave,  or  light  between  a  transmitter 
and  receiver  is  altered,  depending  upon  the  differences 
in  the  respective  transport  properties  between  interior 
heterogeneities  and  their  surroundings.  Consequently, 
information  regarding  interior  heterogeneities  is  con¬ 
tained  in  the  signals  detected  at  the  periphery.  Tomo¬ 
graphic  reconstruction  usually  consists  of  linearizing  the 
relationship  between  received  signals  and  the  spatial 
distribution  of  transport  properties  to  mathematically 
recover  the  spatial  distribution  of  properties  which  yields 
the  best  modeled  match  to  the  receiver  measurements. 
The  result  is  a  tomographic  image  of  the  interior  distri¬ 
bution  of  transport  properties  detailing  the  presence  of 
heterogeneities. 

The  feasibility  of  optical  tomography  has  been  best 
demonstrated  using  single-pixel  measurements  of  fre¬ 
quency-domain  photon  migration  (O'Leary  et  ah,  1995; 
Pogue  et  al. ,  1995).  These  point  measurements  consist 
of  launching  sinusoidally  modulated  light  into  a  highly 
scattering  medium,  such  as  tissue  or  a  particulate 
suspension,  and  measuring  the  modulation  phase  delay, 
6 ,  and  amplitude,  M,  of  the  emitted  light  relative  to  the 
incident  source  (Figure  1).  Solution  of  the  inverse¬ 
imaging  algorithm  involves  iteratively  updating  a  map 
of  optical  properties  until  these  “single-pixel”  measure¬ 
ments  of  phase  shift,  6 ,  and  amplitude  demodulation,  M, 
match  model  predictions.  The  propagation  of  intensity 
modulated  light  can  be  predicted  from  the  diffusion 
equation  (Patterson  et  al.,  1989): 
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Figure  1.  Technique  for  making  frequency-domain  measure¬ 
ments  is  depicted  schematically,  overlaid  on  a  logarithmic  plot 
of  photon  fluence  or  density  produced  by  numerical  simulation. 
White  represents  high  photon  density.  The  medium  is  assumed 
to  have  the  typical  scatter  parameters  of  breast  adipose  tissue 
with  two  heterogeneities  that  have  a  10  x  increase  in  absorption. 

-V*Z)(?)  V4>(r,a;)  +  ~y<E>0 r,co)  +  ^(7)  <3>(r,<y)  = 

S(?,cd)  (1) 

where  <£(m>)  is  the  ac  fluence  rate  of  light  at  position  7 
and  modulated  at  frequency,  co;  //a  is  the  absorption 
coefficient  responsible  for  elimination  of  propagating  light 
by  the  process  of  absorption;  c  is  the  speed  of  light  in  the 
medium;  and  D  is  the  optical  diffusion  coefficient,  given 
by 

3[fia(r)  +  fis'(r')] 


(2) 


The  coefficient,  fis',  is  the  isotropic  scattering  coefficient. 
By  employing  the  proper  boundary  conditions  (Haskell 
et  al,  1994;  Patterson  et  al ,  1989),  values  of  the  complex 
ac  fluence,  <b(7vw),  can  be  calculated  numerically  for 
arbitrary  spatial  distributions  of  the  absorption,  jua,  and 
isotropic  scattering,  ft/,  coefficients.  The  modulation 
phase  shift  and  amplitude  at  any  position  can  be  obtained 
directly  from  the  complex  ac  fluence,  <£(?yw)  =  Meie. 

The  diffusion  approximation  describing  the  propaga¬ 
tion  of  light  is  valid  when  photons  are  predominantly 
scattered  rather  than  absorbed  (i.e.,  fia  « This  is 
indeed  the  case  in  most  tissues,  and  the  diffusion 
equation  has  been  successfully  applied  to  many  problems 
involving  light  propagation  in  biological  media.  H owever, 
in  some  cases,  the  more  general  theory  of  radiative 
transfer  must  be  applied  (Ishamaru,  1978).  In  addition, 
because  the  absorption  coefficient  is  related  to  the 
product  of  extinction  coefficients  and  the  concentration 
of  light  absorbing  compounds  (i.e.,  fia  =  e[C\),  optical 
tomography  also  gives  added  local  information  on  (bio-) 
chemical  composition.  In  the  case  of  tissues,  the  domi¬ 
nant  light  absorbing  species  are  melanin,  water,  and 
deoxy-  and  oxyhemoglobin. 

Typically,  optical  tomography  using  frequency-domain 
photon  migration  measurements  has  been  used  to  recon¬ 
struct  tomographic  images  in  models  or  phantoms  which 
mimic  the  scattering  properties  of  tissues.  O'Leary  and 
co-workers  (1995)  have  successfully  employed  the  first- 
order  Bom  scattering  approximation  to  linearize  the 
relationship  between  detected  fluence  and  spatial  distri¬ 
butions  of  absorption,  fia,  and  isotropic  scattering,  jus\ 
coefficients  to  reconstruct  images  from  experimental, 
near-infrared  measurements.  The  linearization  is  valid 
for  small  changes  in  optical  properties  at  vanishingly 
small  point  volumes.  Pogue  and  co-workers  (1995)  have 
relaxed  these  restrictions  and  have  also  successfully 
reconstructed  images  from  experimental,  frequency- 
domain  photon  migration  measurements.  Their  ap¬ 
proach  uses  an  iterative  approach  based  upon  Newton's 
method  to  converge  upon  maps  of  absorption  and  scat¬ 
tering  by  minimizing  the  differences  between  model 
predictions  and  measured  fluences.  These  studies  pro¬ 
vide  promising  results  for  the  development  of  optical 
tomography  for  biomedical  imaging  (Jiang  et  al. ,  1996; 
Paulsen  and  Jiang,  1996).  Although  future  improve¬ 
ments  to  the  inverse  imaging  solution  are  certain,  two 
stumbling  blocks  remain:  (i)  instrumentation  for  rapid 
registration  of  ac  fluence  signals  and  (ii)  the  availability 
of  endogenous  contrast  for  biomedical  optical  tomogra¬ 
phy. 

Measurements  of  ac  Fluence  for  Tomographic 
Reconstructions.  Frequency-domain  measurements  of 
photon  migration  have  been  demonstrated  to  frequencies 
as  high  as  the  gigahertz  range  (Madsen  et  al .,  1994; 
Fishkin  et  al .,  1996,  1997).  However,  the  attenuation  of 
the  modulation  envelope  increases  rapidly  at  these  higher 
frequencies  due  to  the  order  of  nanosecond  photon  “time- 
of-flight”  in  tissues  (Chance  et  al. ,  1988).  This  increased 
attenuation  coupled  with  decreased  detector  performance 
at  higher  frequencies  limits  the  range  of  useful  modula¬ 
tion  frequencies  for  most  biomedical  optical  imaging 
applications  to  between  30  and  200  MHz  (Sevick  et  al., 
1992).  Typically,  transmitters  consist  of  laser  diodes 
biased  above  the  lasing  threshold  and  modulated  at 
frequency,  f,  by  the  addition  of  an  rf  signal  (Chance  et 
al,  1990;  Thompson  et  al,  1992).  Alternatively,  constant 
intensity  sources  can  be  used  when  coupled  with  an 
acousto-optic  modulator  or  Pockels  cell  assembly  (Gratton 


and  Limkeman,  1983;  Piston  et  al,  1989).  Detection  of 
the  modulation  amplitude  and  phase  at  these  high 
modulation  frequencies  can  be  challenging.  Typically, 
photomultiplier  tube  detectors  with  high  gain  are  em¬ 
ployed  in  heterodyne  mode,  whereby  the  voltage  across 
the  electron  multiplier  chain  is  gain  modulated  at  the 
modulation  frequency,  f,  of  the  incident  light  source 
plus  an  offset  frequency,  A/*,  of  100  Hz  or  more.  This 
mixing  results  in  a  current  at  the  anode  that  is  modu¬ 
lated  at  the  frequency,  A f,  but  contains  the  information 
of  modulation  amplitude  and  phase  shift  of  the  optical 
signal  modulated  at  frequency  f  (Spencer  and  Weber, 
1969;  Graton  and  Limkeman,  1983).  Other  possibilities 
in  which  the  mixing  of  the  electronic  signal  occurs  after 
the  photomultiplier  tube,  rather  than  in  its  electron 
multiplier  chain,  are  also  possible.  The  advantage  of 
this  heterodyned  detection  is  that  simple  data  acquisition 
can  be  used  to  acquire  and  evaluate  the  modulation 
amplitude  and  phase  shift  of  the  ac  fluence  detected  at  a 
single  location.  Typically,  measurement  error  in  modula¬ 
tion  phase  and  amplitude  are  0.1°  and  1%,  respectively 
(Barbieri  et  al,  1989,  1990).  The  disadvantage  for 
these  single-pixel  detection  systems  is  that,  for  a  tomo¬ 
graphic  imaging  system,  multiple  detectors  need  to  be 
replicated  or  multiplexed  to  provide  the  number  of  inputs 
required  to  solve  the  inverse  imaging  problem.  Boas 
and  co-workers  (1997)  report  that  often  the  limiting 
error  in  inverse  solutions  for  optical  tomography  arises 
from  uncertainty  in  detector  position  typically  deter¬ 
mined  from  the  positioning  of  the  fiber  optic  coupling 
leading  from  the  detector  to  the  air-phantom  inter¬ 
face.  In  this  paper  we  demonstrate  multipixel  measure¬ 
ments  of  ac  fluence  directed  for  the  development  of 
biomedical  optical  tomography  that  do  not  require  the 
positioning  of  multiple  detectors  or  single  detector  scan¬ 
ning. 

Available  Contrast  for  Tomographic  Reconstruc¬ 
tions  in  Tissues.  Biomedical  optical  tomography  is 
successful  only  if  the  contrast  in  the  optical  properties 
between  heterogeneities  and  their  surroundings  are 
significant  enough  to  cause  alteration  in  measured 
signals.  However,  in  the  work  by  Troy  et  al  (1996),  the 
endogenous  contrast  available  by  tissue  scattering  was 
not  found  to  be  consistently  large  enough  to  alter  the 
propagation  of  light.  While  their  in  vitro  measurements 
were  unable  to  assess  the  contrast  owing  to  absorption 
due  to  the  increased  hemoglobin  associated  with  angio¬ 
genesis,  single-pixel  ac  measurements  suggest  that  the 
degree  of  contrast  offered  by  the  presence  of  a  perfect 
absorber  (pia  — *  <*=)  in  an  unrealistically  absorption-free 
surroundings  (jua  =  0)  may  still  be  too  small  (Sevick- 
Muraca  et  al,  1997).  Measured  phase  shifts  are  altered 
by  less  than  6°  at  200  MHz  and  amplitude  demodulated 
by  a  few  percent  due  to  the  presence  of  a  perfect  light 
absorbing  heterogeneity  located  in  a  lossless,  highly 
scattering  medium.  This  is  the  maximum  possible 
contrast  owing  to  absorption  and  is  physiologically 
unrealistic.  Consequently,  another  mechanism  of  con¬ 
trast  imposed  by  artificial  or  exogenous  methods  may  be 
necessary  (Sevick-Muraca  et  al,  1997). 

Recently,  we  have  demonstrated  that  the  ac  contrast 
offered  by  fluorescence  exceeds  that  offered  by  absorption 
(Sevick-Muraca  et  al ,  1997).  The  additional  mechanism 
for  contrast  arises  from  the  kinetics  of  the  fluorescence 
decay  process.  Once  a  fluorophore  absorbs  excitation 
light  and  undergoes  an  electronic  transition,  it  remains 
in  an  activated  state  for  a  mean  period  of  time,  known 
as  the  lifetime,  before  relaxation  to  the  ground  state  and 
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Figure  2.  Apparatus  for  multipixel  frequency-domain  photon 
migration  imaging. 

re-emission  of  a  fluorescent  photon  occurs.  Typically 
lifetimes,  r,  are  on  the  order  of  nanoseconds,  inducing 
decades  of  fluorescent  phase  shift  change  and  tens  of 
percent  of  fluorescent  amplitude  demodulation  relative 
to  the  incident  light.  The  phase  shift  and  amplitude 
modulation  of  the  propagated  fluorescent  light  relative 
to  the  incident  excitation  source  can  be  predicted  from 
the  coupled  diffusion  equations  for  both  excitation  and 
emission  (fluorescent)  light  propagation  (Patterson  et  aly 
1994;  Sevick-Muraca  and  Burch,  1994;  Hutchinson  et  at ., 
1996;  Reynolds  et  al .,  1997): 


- V*Dx(r)  V$x0»  +  —<£(/»  + 

Cx 

jx ax(r)  ®x(r,co)  =  Sx(r,(o)  (3) 
— ' V-DJr)  V4>m(r,w)  +  m(7»  + 

vff)  3>m(7»  =  #axr7T7^  (4) 

1  +  (c or ) 

where  the  subscript  x  denotes  the  properties  and  fluence 
for  excitation  light  and  subscript  m  denotes  that  for 
emission  light.  The  absorption  coefficient  for  excitation 
light,  //ax,  is  comprised  of  absorption  due  to  the  fluoro- 
phore,  jUaxp  as  well  as  that  owing  to  chromophores  which 
do  not  result  in  emission  of  light.  The  quantum  efficiency 
of  the  fluorophore  is  denoted  by  <p.  The  term  on  the  right- 
hand  side  of  eq  4  constitutes  the  source  term  of  emission 
light  and  assumes  first-order  fluorescent  decay  kinetics. 

Tomographic  inversion  of  fluorescent  ac  measurements 
using  eqs  3  and  4  constitutes  a  more  complicated  problem 
than  the  nonfluorescent  measurements  in  eq  1.  However, 
the  disadvantage  of  added  complexity  may  be  offset  by 
the  increased  contrast  available  if  enough  ac  signal  can 
be  detected  for  phase  shift  and  amplitude  demodulation 
measurements.  In  addition,  since  the  electronic  relax¬ 
ation  from  activated  to  ground  state  can  be  governed  by 
the  local  biochemical  environment,  the  local  value  of 
lifetime  can  provide  spectroscopic  information.  Despite 
the  added  complexity  of  fluorescence  photon  migration, 
two  groups  have  demonstrated  the  feasibility  of  recon¬ 
structing  lifetime,  r,  and  the  product  of  quantum  yield 
and  fluorophore  absorption  coefficient,  0//^  from  simu¬ 
lated  measurements  (O’Leary  et  al. ,  1996;  Sevick-Muraca 
et  al,  1996;  Paithankar  et  al ,  1997).  Experimental 
measurements  remain  yet  to  be  performed  to  confirm  the 
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Figure  3.  Detailed  circuit  schematic  for  the  gain-modulated  image  intensifier. 
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Figure  4.  Graphical  depiction  of  homodyne  multipixel  data  acquisition. 
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feasibility  of  fluorescence  lifetime  imaging  for  biomedical 
optical  imaging  with  fluorescent  contrast  agents.  Agents 
such  as  indocyanine  green  (ICG)  are  already  approved 
by  the  FDA  for  systemic  introduction  into  the  human 
body  and  have  been  shown  to  provide  measurable  signals 
in  tissue  phantom  studies  in  our  laboratory  (Sevick  et 
aly  1997).  In  this  work,  we  extend  multipixel  measure¬ 
ments  to  fluorescence  measurements  with  the  addition 
of  an  interference  filter  to  select  the  ac  component  of  the 
fluorescent  light.  The  goal  of  our  multipixel  fluorescence 
measurements  is  to  demonstrate  enhanced  contrast 
owing  to  fluorescence  over  absorption,  even  with  the 
reduced  fluorescence  signal  level. 

In  the  following  sections,  we  describe  the  instrumenta¬ 
tion  associated  with  multipixel  ac  fluence  measurements 
and  present  ac  measurements  of  diffusing  photon  density 
waves  propagating  in  tissue-like  scattering  media  with 
absorbing  and  fluorescing  heterogeneities.  By  comparing 
measurements  in  the  presence  and  absence  of  these 
heterogeneities,  we  demonstrate  the  contrast  provided 
by  each  mechanism  as  detected  with  multipixel  measure¬ 
ments. 

Experimental  Technique 

Apparatus.  The  experimental  apparatus  for  multip¬ 
ixel  frequency-domain  photon  migration  measurements 
is  shown  schematically  in  Figure  2  (Troy  et  aL,  1997).  A 
collimated  and  circularized  laser  diode  (Melles  Griot 
model  06DLS403,  789  nm,  30  mW;  model  06DLS503,  836 
nm,  40  mW)  illuminates  the  scattering  phantom  with  a 
sinusoidally  amplitude  modulated  signal.  The  diode  is 
maintained  at  a  constant  dc  bias  midpoint  above  thresh¬ 
old  by  a  Melles  Griot  model  06DLD203  laser  diode  driver 
which  also  maintains  the  diode  at  a  constant  tempera¬ 
ture.  A  13  dBm  (20  mW)  sinusoidal  rf  signal  (Marconi 
2022D  frequency  synthesizer)  is  superimposed  on  the 
constant  dc  current  through  a  bias  network  internal  to 
the  laser  diode  (Thompson,  1992),  The  scattered  signal 
from  the  phantom  is  focused  onto  the  photocathode  of 
the  image  intensifier  by  a  50  mm  Nikorr  camera  lens. 
This  lens  is  positioned  so  that  the  18  mm  diameter  of 
the  intensifier  photocathode  has  a  field  of  view  of  ap¬ 
proximately  80  mm.  An  optical  interference  filter  can 
be  inserted  at  the  input  of  the  image  intensifier  to 
separate  fluorescent  and  excitation  light.  The  weak  near- 
infrared  signal  at  the  photocathode  is  converted  to  an 
electronic  image  which  is  amplified  by  the  image  inten- 
sifier’s  microchannel  plate  and  converted  to  a  visible 
optical  image  when  the  amplified  electrons  strike  the 
phosphor  output  plate.  The  18  mm  diameter  optical 
output  of  the  intensifier  is  then  coupled  to  a  V2  in.  CCD 
array  with  a  105  mm  Nikorr  microlens.  Coupling  with 


Figure  5.  Schematic  illustrating  the  20  cm  cubic  tank  used 
for  a  tissue  phantom  and  the  side  illumination  geometry  used 
for  the  presented  measurements. 


a  tapered  fiber  optic  face  plate  would  be  more  efficient 
and  compact,  but  lens  coupling  allows  for  both  the  use 
of  a  mechanical  shutter  in  front  of  the  CCD  array  and 
direct  viewing  of  the  intensifier  output  with  the  naked 
eye  during  alignment  and  testing.  The  CCD  camera 
(Photometries  AT200)  has  a  resolution  of  512  x  512  pixels 
with  a  16  bit  dynamic  range. 

The  heart  of  this  instrument  is  a  gain-modulated  image 
intensifier.  The  intensifier  and  a  detailed  schematic  of 
its  bias  network  are  shown  in  Figure  3.  The  image 
intensifier  is  an  18  mm  diameter  Generation  III  device 
(Litton  Electronic  Devices  model  510-5916-310)  and  is 
similar  to  those  found  in  common  military  and  com¬ 
mercial  night  vision  systems.  Generation  III  image 
intensifiers  are  characterized  by  their  extended  near- 
infrared  response,  high  gain,  and  high  resolution.  This 
particular  image  intensifier  has  a  photocathode  sensitiv¬ 
ity  of  209  mA/W  at  830  nm  and  103  mA/W  at  880  nm,  a 
current  gain  of  27  500,  and  a  resolution  of  45  lp/mm.  An 
adjustable  high-voltage  power  supply  (GBS  Micro  Power 
Supply  model  PS20060500)  provides  all  the  dc  high 
voltages  to  the  image  intensifier  as  well  as  a  phosphor 
screen  current  limiting  protection  circuit  which  shuts 
down  the  intensifier  when  the  incident  light  level  reaches 
a  potentially  damaging  level.  A  sinusoidal  rf  modulation 
voltage  is  superimposed  on  the  normal  dc  high-voltage 
bias  on  the  image  intensifier’s  photocathode  through 
high-voltage  isolation  capacitors.  The  rf  modulation  is 
typically  a  50  V  peak-to-peak  sinusoid  which  is  provided 
by  an  ENI  model  604L  4  W  linear  rf  power  amplifier. 
Because  the  ENI  amplifier  will  drive  an  open  circuit,  the 
50  Q  power  resistor  can  be  removed  from  the  circuit  to 
increase  the  maximum  rf  voltage  which  can  be  supplied 
to  the  intensifier.  Since  the  gain  of  the  image  intensifier 
is  essentially  proportional  to  applied  voltage,  the  result¬ 
ing  gain  of  the  image  intensifier  is  also  modulated  at  the 


Figure  6.  Multi-pixel  image  of  a  homogeneous  phantom  with  836  nm  side  illumination  (as  in  Figure  5)  and  a  100  MHz  modulation 
frequency.  Image  a  is  the  average  or  dc  pixel  intensities;  image  b  is  the  ac  pixel  intensities;  image  c  is  the  modulation  phase  at  each 
pixel;  and  image  d  is  the  modulation  ratio. 


applied  rf.  Given  a  modulated  optical  signal  at  its  input, 
the  image  intensifier  will  act  as  a  multichannel  rf  mixer 
across  its  field  of  view.  Similar  instrumentation  has  been 
developed  for  fluorescence  lifetime  imaging  in  nonscat¬ 
tering  media  by  Lakowicz  and  Berndt  (1991). 

Data  Acquisition  and  Signal  Processing.  In  the 
application  of  frequency-domain  imaging,  the  gain- 
modulated  intensifier  can  be  used  in  one  of  two  modes. 
In  the  heterodyne  mode,  the  image  intensifier  is  modu¬ 
lated  at  a  frequency  f  +  A fy  offset  from  the  laser  source 
by  A f.  The  difference  frequency  must  be  within  the 
frequency  response  of  the  phosphor  output  of  the  intensi¬ 
fier  (~30  Hz)  and  the  readout  rate  of  the  CCD  array.  The 
rf  amplitude  and  phase  at  a  given  pixel  is  then  regressed 
from  the  magnitude  and  phase  of  the  variation  of  the 
pixel  intensity  at  the  frequency  A /.  In  the  homodyne 
mode,  the  intensifier  is  driven  at  the  same  frequency  as 
the  optical  source.  Here  the  image  produced  is  a  steady- 
state  image  with  no  time  variation.  Rf  amplitude  and 
phase  are  determined  by  varying  the  phase  of  the 
intensifier  relative  to  the  phase  of  the  excitation  light 
and  then  fitting  a  sinusoid  to  the  measured  data.  Ho- 


modyning  has  the  advantage  that  the  image  at  the  output 
of  the  intensifier  is  steady-state  so  that  the  read  out  rate 
of  the  CCD  array  can  be  arbitrarily  slow.  For  this  reason 
we  use  the  homodyne  method  when  taking  data.  How¬ 
ever,  since  we  can  see  the  output  of  the  intensifier 
directly,  we  use  the  heterodyne  method  with  A/^  3  Hz 
to  optimize  the  ac  and  dc  drive  voltages  to  the  intensifier 
by  maximizing  the  A f  flicker  at  the  output  of  the 
intensifier. 

In  our  homodyne  technique,  the  phase  delay,  is 
introduced  by  a  frequency  synthesizer  (Programmed  Test 
Sources  model  PTS  310)  which  is  phase  locked  to  second 
synthesizer  (Marconi  model  2022A)  via  a  10  MHz  phase 
reference  signal,  ensuring  that  the  two  synthesizers  will 
operate  at  exactly  the  same  frequency  with  a  constant 
phase  difference.  The  phase  of  the  PTS  310  synthesizer 
relative  to  the  Marconi  can  then  be  controlled  digitally 
from  0  to  360°  via  an  IEEE-488  GPIB  instrument  control 
bus.  A  homodyned  rf  amplitude  and  phase  picture  is 
then  obtained  by  stepping  the  phase  delay  at  a  regular 
interval  from  0  to  360°  and  acquiring  a  CCD  image  of 
the  intensifier  phosphor  screen  at  each  phase  delay,  as 
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Figure  7.  Multi-pixel  image  of  a  phantom  with  three  absorbing 
(as  in  Figure  5)  and  100  MHz  modulation  frequency.  Image  a  is  th 
image  c  is  the  modulation  phase  at  each  pixel;  and  image  d  is  th 
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eads  of  diameters  4,  6,  and  8  mm  with  836  nm  side  illumination 
average  or  dc  pixel  intensities;  image  b  is  the  ac  pixel  intensities; 
modulation  ratio. 


illustrated  in  Figure  4.  A  single  CCD  image  for  a  given 
relative  phase  shift  is  represented  by  a  matrix  with  each 
i/th  element  representing  a  single  pixel.  The  intensity 
at  a  given  pixel  will  vary  over  a  single  sinusoid  cycle  as 
the  phase  delay  is  varied  over  360°  as  can  be  seen  in 
Figure  4.  Once  this  series  of  images  is  recorded,  the 
measurement  is  complete.  A  simple  and  fast  data 
analysis  algorithm  is  then  applied  to  extract  the  pixel 
by  pixel  intensity  and  phase  data.  The  dc  value  is 
computed  by  taking  the  average  pixel  intensity  for  all 
phase  delays.  The  ac  component  of  the  pixel  intensity  is 
then  isolated  by  subtracting  this  dc  value  from  the  total 
pixel  intensity  at  each  phase  delay.  This  pixel  ac 
intensity  as  a  function  of  phase  delay  is  Fourier  trans¬ 
formed  using  a  Fast  Fourier  Transform  (FFT).  Since 
each  pixel  intensity  varies  over  exactly  one  period,  the 
maximum  of  the  digital  frequency  spectrum  calculated 
by  the  FFT  will  be  found  at  the  first  element  of  the  digital 
frequency  spectrum  and  will  correspond  to  the  frequency 
component  of  the  optical  signal  at  the  modulation  fre¬ 
quency  for  a  given  pixel.  This  complex  number  is  then 
used  to  calculate  the  pixel  phase  (6)  and  rf  amplitude 


(A)  by  the  simple  relationships: 


(5) 


(6) 


where  1(f)  is  the  Fourier  transform  of  the  measured 
intensity  versus  phase  delay,  1(0^),  and  N  is  the  number 
of  phase  delays.  Modulation  ratio  is  simply  the  ratio  of 
the  ac  pixel  intensity  to  the  dc  pixel  intensity,  M  —  ac/ 
dc.  The  FFT  algorithm  is  an  efficient  way  to  find  the  ac 
component  of  the  pixel  intensity.  Ac  magnitude  and 
phase  values  can  be  calculated  approximately  100  times 
faster  than  by  using  a  nonlinear  least-squares  fit  of  the 
data  to  a  sinusoid.  The  FFT  algorithm  is  most  efficient 
when  dealing  with  N  =  2n  data  points,  where  n  is  an 
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Figure  8.  Multi-pixel  ratio  images  (presence/absence)  for  three  absorbing  beads  of  diameters  4,  6,  and  8  mm  with  836  nm  side 
illumination  (as  in  Figure  5)  and  a  100  MHz  modulation  frequency.  Image  a  is  the  average  or  dc  pixel  ratios  (dcpresence/deabsence); 
image  b  is  the  ac  pixel  ratios  (acPresence/acabsence);  image  c  is  the  modulation  phase  difference  (presence  —  ^absence)  pixel;  and  image  d  is 
the  modulation  ratio  (MpreSence/Mabsence)* 


integer.  For  that  reason  we  use  either  32  or  64  phase 
steps  when  taking  our  data. 

The  entire  measurement  procedure  is  automated.  A 
personal  computer  steps  the  phase  of  the  PTS  310  via 
IEEE -488  bus  and  acquires  a  CCD  image  at  each  phase 
delay,  64.  Typically,  three  images  are  taken  at  each 
phase  delay  and  averaged.  The  data  are  then  stored  to 
disk  before  stepping  the  phase  delay  and  repeating.  For 
the  measurements  recorded  here,  the  CCD  array  is 
binned  by  a  factor  of  8  for  a  pixel  resolution  of  128  x 
128.  Typical  exposures  for  our  measurements  range  from 
10  to  1000  ms,  but  measurement  times  are  limited  by 
the  array  readout  time.  Newer  CCD  arrays  are  capable 
of  readout  times  on  the  order  of  0.1  s  with  comparable 
resolution  and  pixel  dynamic  range.  The  data  acquisition 
time  for  a  series  of  100  ms  exposures  would  be  ap¬ 
proximately  equal  to  the  (exposure  time  +  frame  transfer 
time)  x  (number  of  phase  delays)  x  (number  of  pictures 
per  phase  delay)  =  (0.1  +  0.1)  x  64  x  3  =  384  s.  This 
corresponds  to  an  equivalent  acquisition  time  for  data 
acquired  using  “single-pixel”  techniques  of  the  order  T 
x  128  x  128,  where  T  is  the  aquisition  time  required  to 


achieve  a  similar  signal  to  noise  ratio  with  the  single¬ 
channel  detector.  Our  laboratory  experience  indicates 
that  at  these  signal  levels,  single-channel  aquisition 
times  would  be  on  the  order  of  10  s/pixel,  excluding 
mechanical  scan  time.  Thus,  the  total  single-pixel  aqui¬ 
sition  time  would  be  on  the  order  of  160  000  s— over  400 
times  slower  than  this  multipixel  acquisition  system. 

Multipixel  Images  and  Discussion 

To  test  the  effectiveness  of  the  multipixel  technique 
for  frequency-domain  photon  migration  imaging,  we 
imaged  a  number  of  tissue  mimicking  scattering  phan¬ 
toms  with  both  absorbing  and  fluorescing  heterogeneities. 
All  multipixel  images  presented  in  this  paper  are  made 
directly  from  128  x  128  pixel  experimental  images  with 
no  smoothing,  interpolating,  or  any  other  image  process¬ 
ing.  The  tissue  phantom  and  side  illumination  geometry 
that  we  used  in  our  experiments  is  shown  in  Figure  5. 
The  phantom  is  a  20  x  20  x  20  cm3  acrylic  tank  which 
was  filled  with  a  0.5%  by  volume  Intralipid  (Kabi 
Pharmacia,  Inc.)  scattering  solution  which  has  an  iso¬ 
tropic  scattering  coefficient  of  approximately  5.7  cm”1  at 
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Figure  9.  Multi-pixel  difference  images  (presence  -  absence)  for  three  absorbing  beads  of  diameters  4,  6,  and  8  mm  with  836  nm 
side  illumination  (as  in  Figure  5)  and  a  100  MHz  modulation  frequency.  Image  a  is  the  average  or  dc  pixel  ratios  (c?cpresence  “  dcabsence); 
image  b  is  the  ac  amplitude  of  the  complex  difference  for  presence  -  absence;  image  c  is  the  modulation  phase  of  the  complex 
difference  for  presence  —  absence;  and  image  d  is  the  modulation  of  the  complex  difference  for  presence  -  absence. 


a  wavelength  of  789  nm  (van  Staveren  et  al. ,  1991).  Used 
clinically  as  an  intravenous  nutrient,  Intralipid  is  a  milky 
looking  fat  emulsion  that  mimics  the  scattering  proper¬ 
ties  of  tissue  (Driver  et  al,  1989;  van  Staveren  et  al ., 
1991;  Flock  et  al,  1992).  Calibration  of  the  spatial  axes 
for  all  of  the  images  is  achieved  by  imaging  a  ruler  on 
the  face  of  the  phantom  before  making  measurements. 

Absorption  Imaging.  Figure  6  illustrates  the  photon 
migration  images  for  a  homogeneous  tissue  phantom  (no 
heterogeneity)  with  836  nm,  40  mW  excitation  and  100 
MHz  modulation.  Figure  6a  shows  the  dc  or  average 
pixel  intensities;  part  b  is  the  ac  amplitude  or  modulation 
envelope  for  each  pixel;  part  c  is  the  pixel  modulation 
phases;  part  d  is  the  modulation  ratio  for  each  pixel  (ac/ 
dc)  or  part  b/part  a.  Moving  from  right  to  left  across  the 
images,  the  dc  and  ac  intensities  as  well  as  modulation 
ratio  all  gradually  decrease  as  the  distance  from  the 
source  on  the  right  side  (at  x  =  8  cm  and  y  =  4  cm)  of  the 
phantom  increases.  Conversely,  the  phase  delay  of  the 
optical  signal  increases  with  increasing  distance  from  the 
source.  These  images  are  consistent  with  a  diffusing 
photon  density  wave  emanating  from  a  point  source.  Also 


note  that  the  data  plots  are  circular  because  the  field  stop 
for  the  multipixel  imager  is  set  by  the  image  intensifier 
which  has  a  circular  cross  section. 

To  quantify  the  noise  level  at  each  pixel,  a  slice  was 
taken  horizontally  across  the  data  for  each  plot  in  Figure 
6.  Since  this  is  a  homogeneous  scatterer,  all  ac  and  dc 
fluence  profiles  should  have  smoothly  varying  spatial 
profiles  which  are  predicted  by  the  diffusion  equation. 
By  subtracting  this  fitted  smooth  profile  from  the  ex¬ 
perimental  data,  one  obtains  a  measure  of  the  noise  at 
each  pixel.  This  noise  can  be  quantified  as  the  standard 
error  by  computing  the  root  mean  square  of  the  difference 
signal  and  dividing  by  the  square  root  of  the  number  of 
pixels  sampled.  Using  this  aproach  on  the  100  MHz  data 
for  the  homogeneous  phantom,  the  standard  errors  are 
±1.5%  for  the  dc,  ±1.8%  for  the  modulation  amplitude, 
±0.2°  for  the  modulation  phase,  and  ±0.6%  for  the 
modulation  ratio. 

The  data  illustrated  in  Figure  7  are  for  a  similar 
experiment,  but  with  three  black  beads  4,  6,  and  8  mm 
in  diameter  suspended  1  cm  from  the  front  (inside) 
surface  of  the  phantom  and  4  cm  from  the  source.  The 
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Figure  10.  Multi-pixel  ratio  images  (presence/absence)  for  three  absorbing  beads  of  diameters  4,  6,  and  8  mm  with  836  nm  side 
illumination  (as  in  Figure  5)  and  a  30  MHz  modulation  frequency.  Image  a  is  the  average  or  dc  pixel  ratios  (dcpresence/dcabsence);  image 
b  is  the  ac  pixel  ratios  (acpresence/ucabsence) ;  image  c  is  the  modulation  phase  difference  (presence  “  ^absence)  pixel;  and  image  d  is  the 
modulation  ratio  (MpreSenceM/absence). 


shadow  from  the  beads  is  faintly  visible  in  all  four 
images.  To  further  illustrate  the  contrast  resulting  from 
the  heterogeneities,  the  bead  images  are  ratioed  to  the 
“absence”  case  illustrated  in  Figure  6  and  presented  in 
Figure  8.  Here  the  modulation  ratio  (MpreSence/Mabsence) 
and  phase  difference  (presence  ~  ^absence)  are  plotted. 
Although  a  “presence”  and  “absence”  case  of  a  disease 
are  not  clinically  feasible,  we  use  this  to  illustrate  the 
contrast  that  would  be  available  for  tomographic  recon¬ 
struction  algorithms  which  would  use  only  the  presence 
data.  The  degree  of  contrast  can  also  be  assessed  by 
taking  the  complex  difference  between  the  “presence”  and 
“absence”  data.  The  difference  signal  can  then  be  thought 
of  as  the  “scattered”  signal  from  the  heterogeneity.  This 
calculation  is  shown  in  Figure  9.  Again,  the  contrast  can 
be  clearly  seen  in  the  signal  scattered  from  the  hetero¬ 
geneity. 

To  show  the  effect  of  changing  modulation  frequency, 
a  presence/absence  image  for  the  same  phantom  is 
depicted  in  Figure  10  with  a  modulation  frequency  of  30 
MHz.  The  effect  of  reduced  modulation  frequency  can 
be  clearly  seen.  As  expected,  the  dc  image  looks  similar 


but  the  ac-dependent  images  are  different.  However,  the 
longer  wavelength  of  the  “photon  density  wave”  associ¬ 
ated  with  the  lower  modulation  frequency  results  in  a 
much  lower  phase  shift  induced  by  the  absorbers.  The 
phase  noise  of  the  detection  process  then  becomes  com¬ 
parable  to  the  phase  shift  caused  by  the  heterogeneities, 
and  the  contrast  may  be  insufficient  for  image  recon¬ 
struction  using  the  presence  data  alone. 

Fluorescence  Imaging.  Heterogeneities  consisting 
of  fluorescent  dyes  were  imaged  in  the  phantom  to 
demonstrate  the  performance  of  the  imager  and  inves¬ 
tigate  the  contrast  offered  by  fluorescence.  The  excitation 
source  for  these  experiments  is  the  789  nm  laser  diode 
modulated  at  60  MHz.  The  emission  light  is  separated 
from  the  excitation  light  using  an  interference  filter  with 
a  10  nm  bandpass  centered  at  830  nm  (CVI  Laser  Corp. 
model  F10-830-4-100).  The  heterogeneities  consist  of 
micromolar  concentrations  of  the  fluorescent  dye  ICG 
contained  in  a  10  mm  tall,  4  mm  diameter  glass  cylinder 
filled  with  0.5%  Intralipid  scattering  solution.  The 
heterogeneity  is  located  1  cm  from  the  front  surface  of 
the  phantom  and  4  cm  from  the  source.  ICG  has  a 


Figure  11.  Multi-pixel  image  of  a  tissue  phantom  with  a  10  mm  tall  x  4  mm  diameter  glass  cylinder  filled  with  a  2  juM  concentration 
of  ICG  in  0.5%  Intralipid.  The  emission  image  is  at  830  nm  with  789  nm  side  illumination  (as  in  Figure  5)  and  a  60  MHz  modulation 
frequency.  Image  a  is  the  average  or  dc  pixel  intensities;  image  b  is  the  ac  pixel  intensities;  image  c  is  the  modulation  phase  at  each 
pixel;  and  image  d  is  the  modulation  ratio. 


reported  absorption  cross  section  at  780  nm  of  e  = 
130  000  (M  cm)"*1,  a  quantum  efficiency  for  780-830  nm 
conversion  of  0  =  0.016,  and  fluorescent  time  decay  of  r 
=  0.56  ns  (Sevick  et  at. ,  1997).  Because  of  the  reduced 
signal  level  associated  with  fluorescence,  we  used  1000 
ms  exposure  times.  Figure  11  shows  the  images  result¬ 
ing  from  a  2  pM.  concentration  of  ICG.  A  2  juM  concen¬ 
tration  of  ICG  is  approximately  50-75  times  below  its 
lethal  level  and  is  on  the  order  of  that  used  clinically  for 
photodynamic  agents.  Because  of  the  inherent  contrast 
offered  by  fluorescence,  the  object  is  clearly  visible 
without  referencing  to  the  absence  case. 

Conclusion 

We  have  demonstrated  an  efficient  method  for  making 
high  resolution  multipixel  measurements  of  frequency- 
domain  diffusing  photon  density  waves  using  a  gain- 
modulated  image  intensifier.  This  high  resolution  is 
critical  in  making  measurements  for  tomographic  repro¬ 
ductions  because  frequency-domain  measurements  are 
essentially  the  near  field  diffraction  pattern  resulting 
from  the  propagation  of  diffusing  photon  density  waves, 


and  just  as  in  near-field  optical  microscopy,  imaging 
resolution  is  directly  related  to  the  measurement  resolu¬ 
tion.  We  obtain  two-dimensional  maps  of  modulation 
magnitude  and  phase  as  well  as  dc  intensity  with  128  x 
128  pixel  resolution  in  under  7  min.  The  equivalent 
single-pixel  data  acquisition  time  would  be  on  the  order 
of  40  h,  not  counting  mechanical  scan  time.  In  addition, 
because  there  is  no  mechanical  scanning,  the  multi¬ 
pixel  measurement  system  has  much  lower  positioning 
errors  which  have  been  identified  as  the  limiting  noise 
factor  in  scanning  single-pixel  measurements  (Boas  et  ah , 
1997). 

The  increased  gain,  resolution,  and  near-infrared 
response  of  Generation  III  image  intensifiers  allow 
detection  and  localization  of  smaller  heterogeneities  at 
greater  depths  in  tissue  mimicking  phantoms.  In  addi¬ 
tion,  the  extended  infrared  response  also  allows  mea¬ 
surements  to  be  made  over  a  wider  range  of  the  near- 
infrared  biological  window.  This  provides  more  spectro¬ 
scopic  data  and  also  allows  for  low  loss  propagation  at 
both  the  excitation  and  emission  wavelength  when  imag¬ 
ing  fluorescence.  In  addition,  any  exogenous  fluorescent 
contrast  agent  operating  in  this  near-infrared  window 


will  not  be  competing  with  endogenous  fluorescence 
which  is  generally  excited  by  ultraviolet  and  blue 
light. 

Currently,  we  are  working  on  adapting  our  excitation 
and  fluorescent  tomographic  inversion  algorithms  to 
accept  the  wealth  of  data  that  we  obtain  from  our 
multipixel  measurements.  This  rapid  data  acquisition 
coupled  with  improved  inversion  algorithms,  ever  in¬ 
creasing  computational  speeds,  and  exogenous  contrast 
agents  will  hopefully  produce  a  new  and  extremely 
powerful  biomedical  imaging  modality. 
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Notation 

r  position  (cm) 

t  time  (s) 

f  modulation  frequency  (1/s  or  Hz) 

A f  offset  in  modulation  frequency  (1/s  or  Hz) 

a)  modulation  angular  frequency  —  2n f  (rad/s) 

c  speed  of  light  (cm/s) 

i  imaginary  number,  (-1)V2 

jua(r)  absorption  coefficient  or  inverse  mean  free  path 
(1/cm) 

/^axC r)  absorption  coefficient  at  excitation  wavelength 

(1/cm) 

/MaxfCr)  fluorophore  absorption  coefficient  at  excitation 
wavelength  (1/cm) 

/*amCr)  absorption  coefficient  at  emission  wavelength  (1/ 
cm) 

/ is(r )  isotropic  scattering  coefficient  or  inverse  mean 

free  path  (1/cm) 

jusx'(r )  isotropic  scattering  coefficient  at  excitation  wave¬ 

length  (1/cm) 

/*sm'(r)  isotropic  scattering  coefficient  at  emission  wave¬ 
length  (1/cm) 

<f>(r ,  t)  fluence  or  angle  integrated  photon  flux  (W/cm2) 

ZXr)  optical  diffusion  constant  (cm) 

Sx(r,  t)  power  density  source  at  excitation  wavelength 
(W/cm3) 

0  quantum  efficiency  for  excitation  to  emission 

conversion 

r  exponential  time  constant  for  fluorescence  (s) 

0  modulation  phase  (deg  or  rad) 

#d  phase  delay  between  intensifier  and  laser  diode 

modulation  (deg  or  rad) 

M  modulation  amplitude  (W/cm2) 

[C]  molar  concentration  (M) 

€  extinction  coefficient  [(M  cm)"1] 
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Abstract 

An  RF  phase  sensitive  camera  is  used  to  acquire  frequency 
domain  photon  migration  images  of  the  fluorescent  contrast 
agent  indocyanine  green  (ICG)  in  diseased  canine  mammary 
tissue  both  in  vivo  and  ex  vivo.  Using  this  system,  images  of 
100  MHz  modulation  amplitude  and  phase  are  obtained  for 
128  x  128  pixels. 

I.  INTRODUCTION 

The  success  of  frequency  domain  biomedical  optical 
imaging  depends  upon  the  ability  to  measure  changes  in 
photon  migration  due  to  differences  in  the  optical  prop¬ 
erties  of  normal  and  diseased  tissues.  In  fact  several  in¬ 
vestigators  have  successfully  demonstrated  reconstructed 
images  based  upon  differences  in  absorption  and  scatter¬ 
ing  for  idealized  tissue  phantoms  [1.2].  However,  the  dif¬ 
ferences  between  normal  and  diseased  tissue  may  not  be 
consistently  large  enough  for  sufficient  detection  [3].  Con¬ 
sequently.  we  have  embarked  on  studies  to  investigate  ex¬ 
ogenous  contrast  agents  using  multi-pixel  frequency  do¬ 
main  measurements.  In  this  study  we  examine  the  fluo¬ 
rescence  from  ICG  in  diseased  canine  mammary  tissue. 

II.  EXPERIMENTAL  PROCEDURE 

Experimental  measurements  are  performed  using  the 
apparatus  shown  in  Figure  1  [4].  The  light  source  is  a 
30  mW,  789  nm  laser  diode  (circular  and  collimated), 
driven  above  its  lasing  threshold  by  a  constant  current 
source  and  amplitude  modulated  by  a  13  dBm  RF  sig¬ 
nal  supplied  by  a  Marconi  frequency  synthesizer.  For  the 
images  presented  here,  we  use  a  100  MHz  modulation 
frequency.  The  RF  modulated  light  is  delivered  to  the 
tissue  using  a  1  mm  optical  fiber.  The  light  is  allowed 
to  spread  in  its  natural  cone  from  the  fiber  to  a  diame¬ 
ter  of  2.5  cm  to  provide  more  uniform  illumination.  The 
light,  propagates  through  the  tissue  in  a  transillumination 
geometry  and  is  collected  by  the  RF  phase  sensitive  cam¬ 
era.  The  image  intensifier  of  the  phase  sensitive  camera 
is  gain  modulated  at  the  same  phase  locked  frequency  as 
the  laser  diode.  This  produces  a  homodyne  image  of  the 
phantom  where  the  intensity  is  related  to  both  the  av¬ 
erage  intensity  as  well  as  the  modulation  amplitude  and 


phase  of  the  light  emitted  by  the  tissue.  By  taking  a  se¬ 
ries  of  images  with  different  intensifier  phases  relative  to 
the  diode  modulation  phase,  the  individual  pixel  average 
intensity,  modulation  amplitude,  and  modulation  phase 
can  be  determined  efficiently  with  a  Fast  Fourier  Trans¬ 
form  algorithm.  This  technique  is  described  in  more  de¬ 
tail  in  Reference  [4].  Images  of  emission  light  can  be 
obtained  by  inserting  an  appropriate  bandpass  filter  in 
front  of  the  input  to  the  image  intensifier. 


FIG.  1.  Schematic  of  the  multi-pixel  frequency  domain  ap¬ 
paratus. 


The  canine  subjects  in  this  study  are  all  adult,  female 
pets  that  are  about  to  undergo  surgery  to  remove  parts 
of  their  mammary  chain  because  of  palpable  breast  nod¬ 
ules.  Many  of  these  lumps  are  suspected  to  be  naturally 
occurring  breast  cancer.  Thus,  this  model  affords  an  ex¬ 
cellent  opportunity  to  study  the  pharmacokinetics  and  in 
vivo  optical  properties  of  a  fluorescent  contrast  agent  in 
naturally  occurring  diseased  mammary  tissue. 

The  clinical  procedure  is  as  follows.  The  subject  is 
prepared  for  surgery,  which  includes  shaving  the  area 
around  the  mammary  glands.  The  dog  is  anesthetized 
and  moved  into  the  surgical  suite.  While  the  surgical 
team  monitors  the  dog’s  health,  the  RF  phase  sensitive 
camera  is  clamped  to  the  area  with  suspect  breast  lumps. 
The  imager  is  shown  clamped  to  a  canine  breast  in  Fig¬ 
ure  2.  Measurements  of  excitation  and  emission  fluence 
are  made  before  and  after  injection  of  a  clinical  dose  of 


ICG.  After  acquiring  the  in  vivo  data,  the  imager  is  re¬ 
moved  from  the  operating  room  and  a  mastectomy  is 
performed.  Further  ex  vivo  images  are  then  acquired  on 
the  freshly  excised  diseased  tissue.  The  tissue  is  then 
frozen  and  sectioned  for  histology. 


FIG.  2.  Photograph  of  the  apparatus  used  to  clamp  canine 
mammary  tissue  for  transillumination  measurements. 


III.  EXPERIMENTAL  RESULTS 

A  representative  set  of  fluorescent  emission  images  of 
the  right  fifth  (most  posterior)  mammary  gland  is  shown 
in  Figure  3.  This  subject  is  a  11  year  old  miniature  poo¬ 
dle  weighing  7.5  kg.  The  total  exposure  time  for  this  128 
x  128  pixel  modulation  image  is  approximately  1  minute. 
This  data  was  acquired  approximately  25  minutes  after 
systemic  injection  of  a  1.6  cc  0.5%  solution  of  pharma¬ 
ceutical  grade  ICG  in  sterile  saline  solution.  The  source 
is  centered  at  the  lower  edge  of  the  frame  on  the  opposite 
side  of  the  tissue  (transillumination  geometry,  thickness 
1.9  cm).  The  upper  edge  of  the  clamped  breast  tissue  is 
marked  by  the  smooth  transition  from  grey  to  black  near 
1.5  cm  on  the  vertical  axis.  The  raised  tissue  at  3  cm 
on  the  horizontal  axis  is  caused  by  the  suture  which  is 
used  to  support  the  tissue  during  clamping.  The  area 
of  low  phase  around  the  suture  in  Figure  3c.  is  believed 
to  result  from  the  fluorescence  of  a  small  area  of  bleed¬ 
ing  around  the  suture.  The  differences  between  the  in 
vivo  average  (DC)  intensity  images,  modulation  ampli¬ 
tude  images,  and  modulation  phase  images  demonstrate 
that  more  information  can  be  obtained  with  frequency 
domain  measurements  than  with  a  simple  measurement 
of  average  fluorescent  intensity. 
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FIG.  3.  In  vivo  image  of  the  emission  light  from  ICG  at 
830  nm  from  the  right  fifth  mammary  gland  of  a  miniature 
poodle.  The  excitation  light  is  789  nm  with  100  MHz  mod¬ 
ulation.  a)  is  the  average  (DC)  intensity,  b)  is  the  AC  or 
modulation  amplitude,  c)  is  the  modulation  phase,  and  d)  is 
the  modulation  ratio  (AC/DC). 


In  order  to  quantify  the  spatial  uptake  of  the  dye  in 
the  breast  tissue  a  series  of  average  intensity  images  were 
taken  at  periodic  intervals  for  several  minutes  after  injec¬ 
tion  of  the  dye.  The  geometry  and  subject  are  identical 
to  Figure  3.  In  Figure  4  we  present  a  time  histogram 
created  by  taking  a  horizontal  slice  through  each  of  the 
time  images  at  around  1  cm  (through  the  area  of  high 
intensity).  The  data  starts  3  minutes  after  injection  and 
runs  for  23  minutes.  The  fluorescent  intensity  can  be 
clearly  seen  rising  quickly  to  a  peak  at  5  minutes  after 
injection  and  then  dropping  slowly  thereafter  as  the  dye 
washes  out. 


Tmt  after  niection  (nmutes) 


FIG.  4.  Histogram  showing  uptake  and  wash  out  of  ICG  in 
a  mammary  gland  in  vivo.  The  sub  figure  in  the  upper  left 
corner  shows  the  physical  location  of  the  displayed  slice. 

A  representative  ex  vivo  frequency  domain  photon  mi¬ 
gration  image  is  shown  in  Figure  5.  The  excised  tissue  is 
from  the  right  fifth  mammary  gland  of  a  9  year  old,  fe¬ 
male  golden  retriever  weighing  38  kg.  This  image  is  taken 
from  the  skin  side  through  2  cm  of  tissue  in  a  transillumi¬ 
nation  geometry  with  the  source  centered  in  the  frame.  A 
palpable  lump  is  present  in  the  center  of  this  image.  The 
shadow  of  the  lump  (approximately  1  cm  in  diameter) 
can  be  seen  in  the  phase  and  modulation  ratio  images  al¬ 
though  it  is  not  visible  in  the  image  of  DC  intensity.  Once 
again  this  demonstrates  the  added  information  obtained 
using  frequency  domain  photon  migration  imaging. 


FIG.  5.  Ex  vivo  image  of  the  789  nm  excitation  light  from 
the  right  fifth  mammary  gland  of  a  golden  retriever  with  100 
MHz  modulation,  a)  is  the  average  (DC)  intensity,  b)  is  the 
AC  or  modulation  amplitude,  c)  is  the  modulation  phase,  and 
d)  is  the  modulation  ratio  (AC/DC). 


Figure  6  shows  a  time  histogram  of  vertical  slices 
through  the  center  of  the  data  shown  in  Figure  5.  The 
shadow  from  the  tumor  can  be  seen  in  the  middle  of  the 
spatial  axis.  However,  a  definite  shift  in  brightness  from 
bottom  to  top  can  be  seen  as  time  progresses.  We  spec¬ 
ulate  that  this  shift  in  average  intensity  is  due  to  blood 
drainage  in  the  tissue  which  is  clamped  vertically  for  this 
measurement. 


FIG.  6.  Time  histogram  showing  change  in  average  inten¬ 
sity  of  an  ex  vivo  tissue  sample.  The  sub  figure  in  the  upper 
left  corner  shows  the  physical  location  of  the  displayed  slice. 
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IV.  CONCLUSION 

Because  of  its  ability  to  make  rapid  and  highly  sensi¬ 
tive  measurements  of  frequency  domain  diffusing  photon 
density  waves,  an  RF  phase  sensitive  camera  is  proving  to 
be  a  useful  tool  for  clinical  measurements.  In  addition, 
both  the  in  vivo  and  ex  vivo  frequency  domain  images 
demonstrate  that  there  is  important  information  present 
in  the  RF  modulation  envelope  that  is  not  present  in  the 
average  intensity  images.  Continuous  improvements  are 
being  made  to  the  apparatus  to  increase  sensitivity  and 
image  quality.  Clinical  images  allow  for  the  study  of  im¬ 
portant  contrast  agent  parameters  such  as  pharmacoki¬ 
netics,  lifetime,  and  specificity.  Also,  improved  clinical 
images  can  be  used  as  inputs  into  reconstruction  algo¬ 
rithms. 
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Abstract 

Using  a  finite  difference  approach  to  solve  the  coupled  dif¬ 
fusion  equations  for  the  excitation  and  emission  light,  tomo¬ 
graphic  maps  of  absorption,  fluorescence  quantum  efficiency, 
and  lifetime  were  reconstructed  from  simulated  detector  mea¬ 
surements  of  diffuse  photon  density  waves  made  at  the  pe¬ 
riphery  of  a  tissue  phantom. 

I.  INTRODUCTION 

With  the  development  of  near-infrared  (NIR)  laser 
diodes,  the  synthesis  of  fluorescent  dyes  with  excitation 
and  emission  spectra  in  the  NIR  wavelength  regime  has 
accelerated  in  the  past  decade  for  microscopy  applica¬ 
tions.  Owing  to  the  window  of  low  absorbance  in  tissues 
in  this  wavelength  regime,  an  opportunity  also  exists  for 
the  deployment  of  fluorescent  dyes  as  in  vivo  diagnostic 
agents. 

The  detection  of  disease  via  the  use  of  fluorescent  dyes 
and  photodynamic  drugs  as  contrast  agents  has  been  pro¬ 
posed  by  several  investigators.  However,  a  long-standing 
problem  has  been  the  low  uptake  into  diseased  tissues 
providing  insufficient  contrast  for  detection  of  disease. 
While  targeted  delivery  may  improve  uptake  ratios  and 
contrast  imparted  by  these  agents,  contrast  agent  tar¬ 
geting  approaches  have  been  rather  elusive  even  in  con¬ 
ventional  imaging  modalities  such  as  magnetic  resonance 
imaging  and  x-ray  computed  tomography.  Investigators 
have  sought  to  improve  optical  contrast  independent  of 
uptake  ratio  by  employing  dyes  which  fluoresce  differ¬ 
ently  in  diseased  and  normal  tissues. 

One  difficulty,  however,  lies  in  using  multiply  scattered 
light  detected  at  the  air-tissue  interface  to  reconstruct  an 
image  which  differentiates  internal  tissues  on  the  basis  of 
uptake  and  fluorescent  properties  of  quantum  efficiency 
and  lifetime.  Frequency  domain  measurements  of  modu¬ 
lation  phase  and  amplitude  were  conducted  at  the  surface 
of  simulated  tissue  phant  oms  in  order  to  demonstrate  the 
capacity  for  fluorescent  quantum  efficiency  and  lifetime 
imaging.  The  modulation  phase  and  amplitude  of  the 
detected  fluorescent  wave  relative  to  the  incident  excita¬ 
tion  source  provides  information  about  fluorescent  quan¬ 
tum  efficiency  and  lifetime.  However,  due  to  the  intense 
multiple  scattering  properties  of  tissues,  contributions  of 


excitation  and  fluorescent  light  propagation  characteris¬ 
tics  are  also  incorporated  into  the  measurements. 

Using  a  mathematical  model  of  diffuse  light  propa¬ 
gation  and  fluorescence  generation  within  random  me¬ 
dia,  the  presence  and  location  of  simulated  diseased  tis¬ 
sue  heterogeneities  within  a  tissue  volume  was  detected. 
These  computations  suggest  that  multiple  source  and  de¬ 
tection  schemes  need  to  be  utilized  for  biomedical  diag¬ 
nostic  imaging  and  spectroscopy  of  fluorescence  quantum 
efficiency  and  lifetime.  The  method  is  similar  to  that 
used  by  Yorkey  et  al  [1]  for  electrical  impedance  tomog¬ 
raphy.  Pogue  et  al.  [2]  for  absorption  and  scattering  re¬ 
constructions,  and  Paithankar  et  al.  [3]  who  performed 
reconstructions  on  both  fluorescent  lifetime,  r,  and  flu¬ 
orescent  yield  (the  product  of  the  absorption  coefficient 
due  to  fluorophores  and  quantum  efficiency.  <ppaxf)-  In 
addition  to  these  types  of  reconstructions,  fluorescent 
lifetime  imaging  has  been  performed  using  alternative 
approaches  [4,5]. 

II.  FORWARD  PROBLEM 

The  forward  problem  generates  simulated  exterior  de¬ 
tector  data.  This  data  is  then  used  by  the  inversion  al¬ 
gorithm  to  reconstruct  optical  property  maps  of  the  inte¬ 
rior.  The  approach  consists  of  a  two  dimensional,  simu¬ 
lated  tissue  phantom  with  a  4  x  4  cm2  surface  (Figure  1). 
The  phantom  is  discretized  into  either  a  17  x  17  or  a 
33  x  33  grid.  An  optical  heterogeneity  is  placed  in  the 
phantom  and  is  either  one  pixel  (for  the  17  x  17  grid)  or 
nine  pixels  (for  the  33  x  33  grid)  in  size. 

The  source  and  the  detector  locations  were  placed  one 
pixel  in  from  the  edge  of  the  phantom  to  overcome  the 
effects  of  the  zero  fluence  boundary  condition.  A  source 
is  centered  on  each  side  of  the  phantom  for  a  total  of  4 
sources.  Our  approach  to  the  forward  problem  consists 
of  solving  the  coupled  diffusion  equations  using  finite  dif¬ 
ferences  in  order  to  obtain  detector  data  at  either  56  (17 
x  17  grid)  or  120  (33  x  33  grid)  locations  around  the 
periphery  of  the  phantom  for  each  of  the  4  sources. 

To  solve  the  forward  problem,  the  optical  properties  of 
the  background  and  the  object  are  input  as  parameters 
into  the  program  and  then  the  fluence  at  each  grid  point 
is  calculated  using  MUDPACK,  a  finite  difference,  multi- 


grid  software  package  [6].  The  fluence  values  at  the  detec¬ 
tor  positions  are  then  used  to  calculate  the  modulation 
phase,  and  the  log  of  the  modulation  amplitude, 

Jac  (superscripts  x  and  ???  represents  the  excitation  and 
emission  light,  respectively). 


0T =  tan"1 


’  I  w4>xm' 
Re$T'™ 


I ac  =  lo6io  . 


(2,1) 

(2.2) 


During  the  simulation,  only  one  source  is  active  at  a 
time.  A  total  of  four  simulated  experiments  constitute 
the  data  to  be  acquired  experimentally.  In  all  the  sim¬ 
ulated  experiments,  the  sources  are  modulated  at  a  fre¬ 
quency  of  100  MHz.  Typical  forward  calculations  on  a 
Ultra  2  Sun  workstation  take  approximately  1.3  s. 

After  detector  data  is  evaluated  for  Qr'm  and  JT/c  ' 
simulated  Gaussian  noise  with  a  standard  deviation  of 
0.1°  for  phase  and  1%  for  modulation  amplitude  is  added 
using  a  MATLAB  routine.  The  Gaussian  noise  is  added 
to  account  for  the  noise  which  would  be  encountered  in 
a  real  experiment.  The  simulated  forward  results  with 
added  noise  are  then  input  to  the  inversion  algorithm. 


Source 

and 

Detector 


FIG.  1.  Schematic  of  the  two-dimensional  simulated  tissue 
phantom  showing  the  location  of  the  4  sources  and  56  detec¬ 
tors  around  the  perimeter  for  a  17  x  17  grid  of  16  cm2  area. 
The  small  square  object  (one  pixel  in  size)  inside  the  phantom 
represents  a  simulated  diseased  tissue  volume. 


III.  INVERSE  PROBLEM 


The  solution  to  the  inverse  problem  requires  recon¬ 
structing  images  of  the  interior  volume  from  exterior 
measurements  of  9T%m  and  IA™  •  The  inversion  algorithm 
is  conducted  in  two  parts  on  the  basis  of  (i)  absorption 
at  the  excitation  wavelength  due  to  fluorophores,  /iax/ , 


and  (ii)  quantum  efficiency,  <j>,  or  lifetime,  r,  at  the  emis¬ 
sion  wavelength.  The  two  parts  are  similar  except  the 
reconstruction  based  on  //ai/  uses  detector  data  at  the 
excitation  wavelength  while  the  other  part  uses  data  at 
the  emission  wavelength. 

Figure  2  shows  the  flow  diagram  for  inversions. 
First,  a  uniform  guess  for  the  optical  property  map  is 
given.  Typical  values  are  those  found  in  normal  tissue. 
The  forward  solution  is  then  calculated  to  obtain  6X  and 
1XAC  detector  data  for  the  excitation  light.  Next,  two 
Jacobian  matrices,  J(0T,^ax/)  and  J {lxAc^*xj)'  are  con“ 
structed  where  each  entry  of  the  matrix  is  defined  as 
jij  =  dfff  /dfi axj ,j  and  jij  =  dVACJd\iaxSj ,  respectively. 
Each  element  of  the  matrix  describes  the  response  of  the 
source-detector  pair  at  position  i  to  changes  in  //ar/  at 
each  pixel  j.  The  partial  derivatives  are  numerically  ap¬ 
proximated  by  calling  the  forward  solution  a  second  time 


FIG.  2.  Flow  diagram  for  /Jax/  inversions  based  on  photon 
migration  measurements  at  the  excitation  wavelength. 


for  a  1%  (17  x  17  grid)  or  a  0.5%  (33  x  33  grid)  in¬ 
crease  in  the  original  pixel  value  of  //flx/  and  subtract¬ 
ing  the  difference  (i.e.  dOT/d/jaxf  «  A6T/A{jaxf  and 

The  values  of  fiarf  are  updated  using  an  iterative 
Newton-Raphson  technique, 
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where  A/7flx/  provides  a  small  change  for  updating  //flx/ . 
The  value  £  is  a  regularization  parameter  which  is  mul¬ 
tiplied  by  the  identity  matrix,  I.  Because  the  Jacobian 
matrices  are  ill  conditioned  due  to  the  small  sensitivity 
of  jjasJ  far  away  from  the  source  and  detector,  the  regu¬ 
larization  parameter  compensates  by  making  the  matri¬ 
ces  more  diagonally  dominant.  The  parameter  £  is  then 
adjusted  by  a  Marquardt-Levenberg  algorithm  at  every 
iteration  [2].  In  order  to  solve  the  linear  algebraic  equa¬ 
tions  (Equation  3.1)  for  Ajjaxf ,  a  LU  decomposition  back 
substitution  method  is  used.  The  updated  optical  prop¬ 
erty  map  is  then  found  by  adding  A/7ax/  to  the  values 
from  the  previous  iteration. 

The  forward  solution  is  re-calculated  using  the  new 
updated  values  in  order  to  determine  the  reconstruction 
error, 
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where  k  is  the  source  number,  n  is  the  total  number  of 
detectors  (56  for  the  17  x  17  grid  or  224  for  the  33  x 
33  grid),  0cra/  and  /^c.  are  the  predicted  detector  data 
and  6*bf  and  lxACobt  are  the  simulated  experimental  de¬ 
tector  data  for  the  excitation  light.  The  values  of  <tq  and 
t tjac  are  the  typical  standard  deviations  of  noise  in  the 
measurements  and  are  taken  to  be  0.1°  and  1%.,  respec¬ 
tively.  Every  node  on  the  grid  will  yield  an  equation  and 
the  only  known  variables  are  at  the  detector  positions. 
Since  there  are  many  more  nodes  than  detectors,  the  in¬ 
version  scheme  is  under-constrained  and  not  guaranteed 
to  converge  on  the  actual  optical  property  map. 

The  entire  procedure  of  iteratively  adjusting  fiaxf  to 
minimize  the  xl  error  continues  until  a  predetermined 
convergence  criterion  is  met  or  a  maximum  of  50  itera¬ 
tions  is  reached.  The  criterion  is  met  if  any  of  the  fol¬ 
lowing  quantities:  (i)  the  value  of  (ii)  the  absolute 
change  in  \£,  or  (iii)  the  relative  change  in  become 
lower  than  1.0  xl0~5.  Typical  times  to  complete  the  in¬ 
verse  solution  are  approximately  45  minutes  for  the  17  x 
17  grid  and  approximately  4  hours  on  the  33  x  33  grid 
on  a  Ultra  2  Sun  workstation. 

If  only  a  map  of  /iflx/  is  desired,  the  program  stops. 
Otherwise,  the  map  of  pfll/  is  used  in  order  to  compute 


a  second  map  of  either  <f>  or  r  from  detector  data  at  the 
emission  wavelength  (Figure  3).  Currently,  the  compu¬ 
tation  for  maps  of  <p  and  r  are  conducted  in  two  sepa¬ 
rate  programs  in  order  to  evaluate  their  individual  per¬ 
formances. 


FIG.  3.  Flow  diagram  for  <j>  or  t  inversions  based  on  photon 
migration  measurements  made  at  the  emission  wavelength. 

The  algorithm  for  the  emission  loop  is  similar  to  the 
one  described  above.  First,  the  forward  solution  is  cal¬ 
culated  to  obtain  0m  and  I%c  measurements  at  the  de¬ 
tector  positions  for  the  emission  light  using  the  homo¬ 
geneous  map  described  above  except  with  the  updated 
values  of  pax/.  Next,  two  Jacobian  matrices  are  con¬ 
structed:  3(6m}<j>)  and  3(I™c,<j>)  for  the  reconstruction  of 
<p,  or  3(6™  ,r)  and  3(I™c,r)  for  reconstruction  of  r.  Each 
element  of  these  Jacobian  matrices  is  defined  as  jij  = 
de^/dtj,  jij  =  dI™Ct jij  -  dOF/dTj,  and  j{j  = 
dl™cjdrj  *  respectively.  Again,  each  element  of  the  Jaco¬ 
bian  matrix  describes  the  response  of  the  source-detector 
pair  at  position  i  to  changes  in  <j>  or  r  at  each  pixel  j.  The 
partial  derivatives  are  approximated  as  described  above 


where  30™ /d<i>  «  A0m/A<£,  dl^Jd^  *  AI%c/A<t>, 
30™  /3r  «  AO™ /At,  and  31™c/dr  «  AI™c/Ar. 

Equations  3.3  and  3.4  provide  updates  of  <f>  and  r  by 
adding  A<£  and  Af,  respectively,  to  the  values  of  <p  and 
r  from  the  previous  iteration, 
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Again,  the  forward  solution  is  re-calculated  using  the  up¬ 
dated  values  in  order  to  determine  the  reconstruction  er¬ 
ror, 
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where  6™al,  I™c  are  the  predicted  detector  data  and 
<a  and  rj/cobf  are  the  simulated  experimental  detector 
data  for  the  emission  light.  The  values  for  a$  and  crjAC 
are  the  same  as  in  the  absorption  reconstruction,  0.1° 
and  1.0%,  respectively.  The  entire  procedure  of  itera¬ 
tively  adjusting  o  or  r  to  minimize  the  \2,  error  contin¬ 
ues  until  the  same  convergence  criterion  as  above  is  met 
or  50  iterations  have  passed. 

Three  different  types  of  simulations  were  separately 
performed.  They  included  reconstructions  on  the  basis 
of  absorption,  i*axr  quantum  efficiency,  6  or  lifetime,  r. 
In  all  three  classes  of  simulated  experiments,  an  optical 
heterogeneity  on  the  17  x  17  grid  was  0.063  cm2  which 
constituted  0.39%  of  the  total  area.  For  the  33  x  33 
grid,  the  heterogeneity  was  0.141  cm2  or  0.88%  of  the 
total  area.  During  the  reconstructions  of  both  <j)  and 
t  unphysically  high  optical  property  values  around  the 
periphery  and  especially  near  the  source  locations  were 
obtained.  These  values  were  replaced  with  the  average 
background  value  with  a  constraint  statement  inside  the 
program. 

For  all  the  reconstructions,  values  of  /jaxJ  range  from 
0.010  to  0.200  cm"1.  These  values  of  absorption  corre¬ 
spond  to  an  ICG  concentration  of  approximately  0.076 
to  1.53  /iM  [7].  These  ICG  concentrations  are  well  below 
lethal  levels  and  are  approximately  5  to  60  times  lower 
than  the  therapeutic  concentrations  currently  adminis¬ 
tered  with  many  photodynamic  agents  [8,9].  The  value 
for  the  isotropic  scattering  coefficient  ,  fj*s ,  was  taken  to 
be  10  cm"3  inside  the  object  and  the  background. 


IV.  RECONSTRUCTION  OF  THE  ABSORPTION 

COEFFICIENT  FROM  EXCITATION  LIGHT 

Figure  4  show  a  reconstruction  based  on  a  10:1  dif¬ 
ference  in  //aj/  between  the  object  and  its  surroundings 
using  the  17  x  17  grid.  As  described  above,  this  ratio  is 
the  current  performance  observed  with  contrast  agents. 
The  values  of  fjQxf  in  the  object  and  background  were 
taken  to  be  0.200  cm”1  and  0.020  cm"3,  respectively. 
Figure  4  shows  the  (a)  ideal  spatial  map,  and  (b)  recon¬ 
structed  image.  Although  the  actual  object  location  is  at 
pixel  (11.7),  the  largest  value  in  the  reconstruction,  0.042 
cm"1,  was  found  at  pixel  (12,7)  which  is  one  pixel  closer 
to  the  detectors  in  the  y-direction.  The  reconstruction  of 
[iaxf  at  location  (11,7)  was  0.035  cm"1  which  is  only  17% 
lower  than  the  pixel  (12,7)  but  75%  larger  than  the  av¬ 
erage  background  value.  Because  the  program  preserves 
the  average  pixel  value,  a  smoothing  of  the  reconstructed 
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FIG.  4.  Reconstructed  spatial  map  for  a  10:1  uptake  ra¬ 
tio  for  an  absorbing  heterogeneity  on  a  17  x  17  grid.  Fig¬ 
ure  (a)  represents  the  actual  spatial  map,  and  (b)  repre¬ 
sents  the  reconstructed  image.  The  actual  values  of  ^Qjc/ 
in  the  object  and  the  surroundings  were  0.200  cm"1  and 
0.020  cm”1,  respectively,  and  the  reconstruction  yielded  val¬ 
ues  of  0.035  cm”1  and  0.020  cm"1,  respectively. 


image  is  observed  causing  the  magnitude  of  fiaxJ  for  both 
the  object  and  the  surroundings  to  be  incorrect. 

Next,  the  performance  of  a  10:1  uptake  ratio  for  a 
reconstruction  using  a  33  x  33  grid  was  investigated 
(Figure  5).  A  more  refined  grid  will  yield  smaller  dis¬ 
cretization  errors  because  the  derivatives  are  better  ap¬ 
proximated.  However,  the  inversion  problem  is  under- 
constrained,  as  discussed  above,  and  is  not  guaranteed 
to  converge  on  a  unique  solution.  Therefore,  a  more  re¬ 
fined  grid  might  not  yield  better  reconstructions  since 
the  number  of  nodes  (equations  to  be  solved)  increases 
quadratically  while  the  number  of  detectors  (known  vari¬ 
ables)  only  increases  linearly. 
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FIG.  5.  Reconstructed  spatial  map  for  an  absorbing  het¬ 
erogeneity  on  a  33  x  33  grid.  Figure  (a)  represents  the  actual 
image,  and  (b)  represents  the  reconstructed  image.  The  ac¬ 
tual  values  of  fjQxJ  in  the  object  and  the  surroundings  were 
0.200  cm-1  and  0.020  cm”1,  respectively,  and  the  reconstruc¬ 
tions  yielded  values  of  0.084  cm”1  and  0.021  cm"1,  respec¬ 
tively. 

The  same  optical  parameters  and  reconstruction  con¬ 
straints  as  with  the  17  x  17  grid  (Figure  4)  were  main¬ 
tained.  Figure  5  shows  (a)  the  actual  spatial  map  and 
(b)  the  reconstructed  image.  This  time  the  reconstruc¬ 
tion  correctly  locates  the  object  since  the  largest  pixel 


value  of  0.084  cm”1  is  found  at  pixel  (25,8).  Although, 
the  magnitude  for  ptaxJ  in  both  the  object  and  the  sur¬ 
roundings  is  incorrect  due  to  the  smoothing  effect,  an 
improvement  over  the  17  x  17  grid  is  seen.  This  recon¬ 
struction  also  converges  after  only  3  iterations  which  is 
faster  than  the  17  x  17  reconstruction  which  took  8  it¬ 
erations. 


V.  RECONSTRUCTION  OF  QUANTUM 
EFFICIENCY  FROM  EMISSION  LIGHT 

The  performance  of  the  inverse  algorithm  for  recon¬ 
structing  maps  of  quantum  efficiency,  <£,  using  detector 
data  at  the  emission  wavelength  is  discussed  in  this  sec¬ 
tion.  For  all  the  inversions  discussed  here,  the  exact  value 
of  inside  the  object  was  input  into  the  initial  homo¬ 
geneous  guess.  This  procedure  causes  the  inversion  al¬ 
gorithm  to  converge  after  one  iteration  in  the  excitation 
loop  and  the  exact  value  of  fiaxS  to  be  used  in  the  emis¬ 
sion  loop  in  order  to  only  investigate  the  reconstruction 
on  <j>.  Also,  the  known  value  of  r  both  inside  the  object 
and  background  were  input  as  part  of  the  initial  guess. 

Figure  6  shows  the  (a)  actual  solution,  and  (b)  recon¬ 
structed  image  using  a  17  x  17  grid.  The  values  of  <j) 
inside  the  object  and  surroundings  were  0.10  and  0.010, 
respectively.  Also,  the  values  of  fiaxf  in  the  object  and 
surroundings  were  0.200  cm”1  and  0.020  cm"1,  respec¬ 
tively.  The  experiment  ended  after  50  iterations  with¬ 
out  meeting  any  of  the  convergence  criteria.  The  optical 
property  map  of  the  reconstruction  shows  that  the  loca¬ 
tion  of  the  object  is  successfully  located  but  the  value  of 
<p  inside  the  object  only  reaches  0.082  (18%  smaller  than 
the  actual  value).  Also,  the  background  pixel  values  are 
noisy  with  an  average  value  of  0.012  which  is  20%  larger 
than  the  actual  value. 

Next,  a  10:1  uptake  of  dye  was  investigated  using  the 
finer  grid  which  is  illustrated  in  Figure  7.  The  same  op¬ 
tical  parameters  as  before  were  used.  The  results  for  the 
reconstruction  are  inconclusive  since  the  optical  property 
map  was  too  noisy  to  obtain  any  useful  information.  The 
reconstruction  yielded  an  average  pixel  value  of  0.051  for 
<p  over  the  entire  map  which  is  approximately  400%  larger 
than  the  actual  background  value.  To  perform  a  recon¬ 
struction  of  this  kind  on  a  33  x  33  grid  more  data  is 
necessary,  either  using  more  detectors  or  more  frequency 
information. 


VI.  RECONSTRUCTION  OF  LIFETIME  FROM 
EMISSION  LIGHT 

The  performance  of  the  inversion  algorithm  for  recon¬ 
structing  lifetime,  r,  using  detector  data  at  the  emission 
wavelength  is  discussed  in  this  section.  Similar  to  the 
reconstructions  of  <j>,  all  the  inversions  in  this  section  use 


the  exact  value  of  fjaxf  inside  both  the  object  and  back¬ 
ground  for  the  initial  optical  property  map. 
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FIG.  6.  Reconstructed  spatial  map  of  fluorescent  quantum 
efficiency,  <£,  on  a  17  x  17  grid  for  a  10:1  uptake  of  dye  inside 
the  heterogeneity.  Figure  (a)  represents  the  actual  image, 
and  (b)  shows  the  corresponding  reconstruction.  The  actual 
values  of  in  the  object  and  the  surroundings  were  0.100  and 
0.010,  respectively,  and  the  reconstructions  yielded  values  of 
0.082  cm  and  0.012  cm-1,  respectively. 

The  lifetime  reconstruction  for  a  10:1  uptake  of  dye, 
was  investigated  on  both  a  17  x  17  and  33  x  33  grid. 
For  all  r  reconstructions,  r  had  an  arbitrary  lifetime  of 
1.00  ns  inside  the  object  and  a  10.0  ns  everywhere  else. 
As  with  <f>  reconstructions,  the  values  of  fiarf  in  the  ob¬ 
ject  and  surroundings  were  0.200  cm-1  and  0.020  cm-1, 
respectively.  The  first  set  of  images  is  for  the  course  grid. 
Figure  8  shows  the  (a)  actual  spatial  map  and  the  (b)  re¬ 
constructed  spatial  map  of  r.  Because  r  inside  the  object 
was  lower  than  in  the  background,  the  object  location  is 
defined  as  the  pixel  with  the  lowest  value  of  r.  Both  the 
location  and  magnitude  of  r  inside  the  object  are  cor¬ 
rectly  found  since  the  reconstructed  value  for  pixel  (11,7) 
is  0.996  ns,  only  0.04%  lower  than  the  actual  value.  The 
reconstruction  for  the  background  converges  to  10.1  ns 
even  though  the  background  contains  a  small  amount  of 


noise  which  is  symmetric  about  the  object.  This  noise  is 
most  likely  an  artifact  from  the  large  fluence  that  emerges 
from  the  four  source  locations. 

The  next  set  of  reconstructions  were  identical  to  the 
last  except  that  a  finer  grid,  33  x  33,  was  used.  Figure  9 
shows  the  (a)  actual  image  along  with  the  (b)  correspond¬ 
ing  reconstruction  for  r.  The  algorithm  does  not  cor¬ 
rectly  reconstruct  the  heterogeneity  but  a  disturbance  in 
the  upper  left  quadrant  where  the  object  is  located  is  ob¬ 
tained.  The  smallest  pixel  value  of  1.00  x  10"5  ns  is 
found  at  grid  number  (26,8)  which  is  one  pixel  location 
above  the  actual  object  location.  Unlike  the  absorption 
case,  the  refined  grid  does  not  yield  better  results,  in  fact, 
this  result  suggests  that  the  algorithm  based  on  lifetime 
cannot  be  used  to  detect  the  exact  location  of  a  hetero¬ 
geneity  containing  a  10:1  uptake  of  dye  using  a  33  x  33 
grid  unless  more  information  is  given. 
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FIG.  7.  Reconstructed  spatial  map  of  fluorescent  quantum 
efficiency,  <j>,  on  a  33  x  33  grid  for  a  10:1  uptake  of  dye  inside 
the  heterogeneity.  Figure  (a)  shows  the  actual  image  and  (b) 
is  the  corresponding  reconstruction.  The  actual  value  of  <t> 
in  the  object  and  background  was  0.100  and  0.010,  respec¬ 
tively.  The  algorithm  was  unable  to  reconstruct  the  object. 
The  average  quantum  efficiency  over  the  entire  map  for  the, 
reconstruction  was  0.053. 


The  final  reconstruction  was  for  a  100:1  uptake  of  dye 
inside  the  object  using  the  finer  grid.  The  values  of  r 
inside  the  object  and  the  backgound  were  kept  the  same 
as  above.  Also,  the  values  of  ytQxJ  in  the  object  was  kept 
the  same,  0.200  cm""1.  This  time  however,  the  value  of 
for  the  the  surroundings  was  set  to  0.002  cm"1. 
Figure  10  shows  (a)  the  actual  image,  and  (b)  the  recon¬ 
struction.  The  smallest  value  of  r,  0.98  ns,  was  found 
at  pixel  location  (25,8).  This  pixel  along  with  its  eight 
surrounding  neighbors  have  an  average  value  of  1.01  ns 
which  is  1.0%  larger  than  the  true  value.  Therefore,  at 
100:1  uptake,  both  the  magnitude  and  the  location  of  the 
object  are  correctly  obtained.  The  same  kind  of  symmet¬ 
ric  noise  about  the  object  that  was  observed  in  the  10:1 
reconstruction  on  the  17  x  17  grid  (Figure  8)  is  observed. 
Again,  the  noise  seems  to  be  a  artifact  of  the  source  lo¬ 
cations. 
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FIG.  8.  Reconstructed  spatial  map  of  fluorescent  lifetime, 
r,  on  a  17  x  17  grid  for  a  10:1  uptake  of  dye  inside  the  het¬ 
erogeneity.  Figure  (a)  shows  the  actual  image  and  (b)  is  the 
corresponding  reconstruction.  The  actual  lifetimes  for  the 
object  and  background  are  1.0  and  10.0  ns,  respectively.  The 
reconstruction  converges  correctly  to  1.0  ns  in  the  object.  The 
background  converges  to  10.067  ns  which  is  0.67%  larger  than 
the  actual  value  of  10  ns. 


FIG.  9.  Reconstructed  spatial  map  of  fluorescent  lifetime, 
t,  on  a  33  x  33  grid  for  a  10:1  uptake  of  dye  inside  the  het¬ 
erogeneity.  Figure  (a)  represents  the  actual  image  and  (b)  is 
the  corresponding  reconstruction. 


VII.  CONCLUSIONS 

Numerical  simulations  were  conducted  in  order  to  ex¬ 
amine  the  resolution  and  accuracy  of  our  reconstruction 
method.  These  results  show  the  potential  of  photon  mi¬ 
gration  for  the  optical  detection  of  a  heterogeneity  inside 
a  scattering  medium.  The  simulated  experiments  based 
on  fiaxJ  can  locate  an  absorbing  heterogeneity  inside  a  tis¬ 
sue  simulating  phantom,  however  the  magnitude  of  fjarf 
inside  the  heterogeneity  is  smaller  than  the  actual  solu¬ 
tion.  The  simulated  experiments  also  show  the  ability  to 
reconstruct  the  interior  optical  property  maps  for  both 
quantum  efficiency  and  lifetime  using  detector  informa¬ 
tion  at  the  emission  wavelength.  Reconstructing  maps  of 
lifetime  not  only  provides  object  location  but  also  may 
provide  information  regarding  the  biochemical  environ¬ 
ment.  This  information  about  the  local  environment  can 
help  to  differentiate  normal  from  diseased  tissue.  Since 
targeted  delivery  of  a  contrast  agent  may  not  always  be 
feasible,  lifetime  sensitive  dyes  which  are  specific  to  dif- 


ferent  environment  al  conditions  found  in  normal  and  dis¬ 
eased  tissue  may  enhance  detection. 

The  noise  for  the  fluorescent  measurements  needs  to 
be  determined  since  it  is  unlikely  to  be  the  same  as  the 
excitation  data  (0.1°  for  6T  and  1%  for  lrAC)'  Also,  both 
inversions  assume  that  the  noise  is  uniform  over  the  en¬ 
tire  region  which  may  not  be  the  actual  situation,  and 
updates  to  the  inversion  scheme,  which  take  the  noise 
as  a  function  of  position  into  account  ,  need  to  be  imple¬ 
mented. 
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FIG.  10.  Reconstructed  spatial  map  of  fluorescent  lifetime, 
r,  on  a  33  x  33  grid  for  a  100:1  uptake  of  dye  inside  the 
heterogeneity.  Figure  (a)  represents  the  actual  image  and 
(b)  is  the  corresponding  reconstruction.  The  reconstructed 
background  has  an  average  value  of  11.1  ns  which  is  11.1% 
from  the  true  value  of  10  ns.  The  reconstructed  lifetime  inside 
the  object  converged  to  an  average  value  of  1.01  ns  which  is 
1.0%.  from  the  true  value  of  1.00  ns. 

To  further  improve  the  reconstructions,  detector  infor¬ 
mation  at  multiple  modulation  frequencies  needs  to  be 
incorporated  into  the  algorithm.  Multiple  frequency  in¬ 
formation  will  provide  more  detector  information  helping 
with  the  under-constrained  nature  of  the  inversion  prob¬ 
lem.  Another  way  to  improve  the  inversion  is  to  solve 
the  forward  problem  on  a  grid  with  a  higher  resolution 


than  the  optical  property  grid.  This  methodology  will 
increase  the  number  equations  while  keeping  the  number 
of  unknowns  constant  thus  making  the  inversion  better 
constrained. 


ACKNOWLEDGMENTS 

This  work  was  supported  in  part  by  the  National 
Institutes  of  Health  (R01CA61413,  R0167176-01  and 
K04CA687374-01). 


[1]  T.  J.  Yorkey,  J.  G.  Webster,  and  W.  J.  Tompkins,  “Com¬ 
paring  reconstruction  algorithms  for  electrical  impedance 
tomography,”  IEEE  Trans.  Biomed.  Eng.  BME-34,  843- 
852  (1987)! 

[2]  B.  W.  Pogue,  M.  S.  Patterson,  H.  Jiang,  and  K.  D. 
Paulsen,  “Initial  assessment  of  a  simple  system  for  fre¬ 
quency  domain  diffuse  optical  tomography,”  Phys.  Med. 
Biol.  40.  1709-1729  (1995). 

[3]  D.  Y.  Paithankar,  A.  U.  Chen,  B.  W.  Pogue,  M.  S.  Pat¬ 
terson,  and  E.  M.  Sevick-Muraca,  “Imaging  of  fluorescent 
yield  and  lifetime  from  multiply  scattered  light  reemitted 
from  random  media,”  Appl.  Opt.  36,  2260-2272  (1997). 

[4]  M.  O'Leary,  D.  Boas,  D.  Li,  B.  Chance,  and  A.  Yodh, 
“Fluorescent  lifetime  imaging  in  turbid  media,”  Optics 
Letters  21,  158-160  (1996). 

[5]  J.  Chang,  H.  Graber,  and  R.  Barbour,  “Luminescence  Op¬ 
tical  Tomography  in  Dense  Scattering  Media,”  JOS  A  A. 
14  (1997). 

[6]  J.  C.  Adams,  “MUDPACK:  Multigrid  portable  FOR¬ 
TRAN  software  for  the  efficient,  solution  of  linear  elliptic 
partial  differential  equations,”  Appl.  Math.  Comput.  34, 
113-146  (1989). 

[7]  E.  Sevick-Muraca,  G.  Lopez,  J.  Reynolds,  T.  Troy,  and  C. 
Hutchinson,  “Fluorescence  and  absorption  contrast  mech¬ 
anisms  for  biomedical  optical  imaging  using  frequency- 
domain  techniques,”  J.  Photochem.  Photobiol.  66,  55-64 
(1997). 

[8]  Photodynamic  Therapy  basic  principles  and  clinical  appli¬ 
cation ,  B.  W.  Henderson  and  T.  J.  Dougherty,  eds.,  (Mar¬ 
cel  Dekker,  Inc.,  1992). 

[9]  S.  L.  Marcus,  “Photodynamic  Therapy  of  Human  Can¬ 
cer,”  IEEE  80,  869-889  (1992). 


Fluorescence  Lifetime  Imaging  and  Spectroscopy 

in  Random  Media 

Tamara  L.  Troy  and  Eva  M.  Sevick-Muraca 

The  Photon  Migration  Laboratory 
School  of  Chemical  Engineering 
Purdue  University 
West  Lafayette,  IN  47907-1283  USA 

http:/ / photon.ecn.purdue.edu/~chepmi  /  ppml.html 


submitted  for  publication  in 

Applied  Fluorescence  in  Chemistry,  Biology  and  Medicine 


TABLE  OF  CONTENTS 


Page 

1.  PART  I  -PRINCIPLES  OF  PHOTON  MIGRATION .  5 

1.1  Time  and  Frequency  Domain  Photon  Migration .  5 

1.2  Transport  Property:  Absorption  Coefficient  .  6 

1.3  Transport  Property:  Scattering  Coefficient .  7 

1.4  Transport  Property:  Anisotropy  Coefficient .  7 

1.5  Governing  Equations  Describing  the  Transport  of  Light .  8 

1.5.1  One-Speed  Approximation  to  the  Transport  Equation .  9 

1.5.2  Diffusion  Equation  Describing  Light  Propagation  in  Tissues  .  9 

1.5.3  Source  Location  and  Boundary  Conditions .  10 

1.6  Non-Radiative  Mechanisms  for  Exogenous  Contrast .  10 

1.6.1  Contrast  Owing  to  Absorption .  11 

1.6.2  Contrast  Owing  to  Dye  with  Fluorescent  Characteristics  ...  12 

1.6.3  Diffusion  Equations  Describing  Excitation  and  Emission  Pho¬ 
ton  Migration .  16 

1.6.4  Potential  Fluorescent  Contrast  Agents  for  Photon  Migration 

Imaging .  17 

2.  PART  II-CONTRAST  FOR  PHOTON  MIGRATION  IMAGING:  FLUO¬ 
RESCENCE  LIFETIME .  19 

2.1  Single  Pixel  Measurements  for  Photon  Migration  Imaging .  19 

2.2  Phantom  Apparatus  .  19 

2.3  Frequency  Domain  Instrumentation  and  Set-Up .  20 

2.3.1  Heterodyne  Detection  Approach  with  Modulated  Laser  Source  20 

2.3.2  Heterodyne  Detection  Approach  with  Pulsed  Laser  Source  .  .  21 

2.3.3  Instrumentation .  21 

2.4  Contrast  Owing  to  Perfect  Absorption  and  Perfect  Uptake  of  Fluores¬ 
cence  .  22 

2.5  Contrast  Owing  to  the  Influence  of  Lifetime .  24 

2.6  Discussion .  24 

3.  PART  III-RECONSTRUCTION  OF  FLUORESCENCE  QUANTUM  EF¬ 
FICIENCY  AND  LIFETIME .  26 


2 


Page 


3.1  INVERSE  IMAGING  PROBLEM .  26 

3.1.1  Forward  Problem .  27 

3.1.2  Solution  to  the  Inverse  Problem .  29 


4.  PART  IV-RECONSTRUCTION  OF  ABSORPTION  AND  FLUORESCENCE 
QUANTUM  EFFICIENCY  AND  LIFETIME  IN  A  SCATTERING  MEDIUM  33 


4.1  Results  for  the  Reconstruction  Algorithm .  33 

4.2  Reconstruction  of  the  Absorption  Coefficient  from  Excitation  Light  .  33 

4.3  Reconstruction  of  Quantum  Efficiency  from  Emission  Light .  34 

4.4  Reconstruction  of  Lifetime  from  Emission  Light .  35 

4.5  Discussion .  36 

4.6  Conclusion .  37 

LIST  OF  REFERENCES .  51 


i 


3 


With  the  development  of  near-infrared  (NIR)  laser  diodes,  the  synthesis  of  fluo¬ 
rescent  dyes  with  excitation  and  emission  spectra  in  the  NIR  wavelength  regime  has 
accelerated  in  the  past  decade  for  microscopy  applications.  Owing  to  the  window  of 
low  absorbance  in  tissues  in  this  wavelength  regime,  an  opportunity  also  exists  for 
the  deployment  of  fluorescent  dyes  as  in  vivo  diagnostic  agents.  Figure  1.1  illustrates 
the  typical  absorbance  spectra  of  tissues  showing  that  at  wavelengths  less  than  650 
nm,  hemoglobin  absorbance  provides  the  predominant  attenuation  of  light  in  tissues, 
and  above  900  nm  water  absorbance  provides  the  predominant  attenuation  of  light  in 
tissues.  In  the  “therapeutic  window”  of  650-900  nm,  a  window  of  low  absorption  ex¬ 
ists  in  which  light  will  be  preferentially  scattered  over  being  absorbed.  Consequently, 
it  is  possible  to  transmit  multiply  scattered  NIR  light  across  several  centimeters  of 
tissue.  In  addition,  it  is  possible  to  excite  NIR  fluorophores  deep  within  tissues  and 
to  collect  the  fluorescence  re-emitted  from  the  air-tissue  interface.  Since  fluorescence 
provides  a  sensitive  means  for  assessing  local  biochemistry  via  changes  in  quantum 
efficiency  and  lifetime,  the  ability  to  diagnose  tissues  based  on  spectroscopic  evalu¬ 
ation  of  lifetime  and  quantum  efficiency  of  exogenous  diagnostic  agents  is  possible. 
Agents  whose  emission  characteristics  vary  with  tissue  pH  [1,  2]  and  pC>2  [3]  have  been 
employed  to  detect  diseased  tissues  by  the  nature  of  differing  fluorescent  properties  as 
well  as  to  provide  diagnostic  information  regarding  the  diseased  tissue  volume.  While 
typical  contrast  agents  employed  for  the  detection  of  diseased  tissues  depend  on  and 
are  limited  by  the  poor  preferential  uptake,  the  alteration  of  fluorescent  properties 
provides  a  unique  mechanism  for  inducing  additional  contrast  [4].  For  example,  using 
time-gating  to  collect  the  long-lived  fluorescence  from  hematoporphyrin  derivative 
(HpD),  Cubbeddu  and  coworkers  [5,  6]  were  able  to  differentiate  normal  tissues  from 
intradermally  or  subcutaneously  implanted  murine  tumors  in  mice.  More  recently,  it 
has  been  reported  that  the  fluorescent  decay  of  HpD  is  appreciably  slower  in  experi¬ 
mental  mice  tumors  than  in  their  surrounding  normal  healthy  tissues.  Consequently, 
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the  use  of  a  fluorescent  dye  may  provide  contrast  owing  to  changes  in  fluorescent 
properties  within  tissue  compartments. 

The  difficulty,  however,  lies  in  understanding  the  use  of  multiply  scattered  exci¬ 
tation  light  to  excite  a  fluorophore  in  the  tissue,  and  secondly,  to  extract  information 
of  lifetime  and  quantum  efficiency  from  the  multiply  scattered  fluorescence  detected 
at  the  air-tissue  interface.  If  the  optical  properties  of  the  tissue  are  spatially  uniform 
and  the  fluorescent  dye  has  constant  fluorescent  properties,  then  the  problem  is  one 
of  fluorescent  lifetime  spectroscopy  in  tissues.  However,  if  the  detection  of  diseased 
tissues  is  to  be  tackled,  the  problem  becomes  one  of  fluorescent  lifetime  spectroscopic 
imaging ,  since  optical  properties  and  fluorescent  properties  of  diagnostic  agents  may 
vary  with  spatial  location. 

In  this  chapter,  frequency  domain  photon  migration  fluorescence  imaging  is  de¬ 
scribed  as  a  method  for  generating  an  optical  map  or  image  of  fluorescent  lifetime  and 
quantum  efficiency  from  exterior  measurements  of  modulation  phase  and  modulation 
amplitude  on  tissues  or  highly  scattering  media.  In  Part  I,  the  theory  behind  the  prop¬ 
agation  of  excitation  light  and  generation  of  fluorescent  light  within  scattering  media 
such  as  tissues  is  presented.  Part  II  describes  experimental  measurements  which  show 
that  the  delay  in  phase  and  decrease  in  amplitude  of  fluorescence  measured  in  sim¬ 
ulated  tissue  phantoms  varies  as  a  function  of  dye  lifetime.  Part  III  describes  the 
general  theory  behind  the  derivation  of  an  optical  property  map  from  measurements 
of  modulation  phase  and  amplitude  of  fluorescent  light  detected  at  the  tissue  surface. 
Part  IV  presents  actual  images  generated  from  simulated  measurements  of  modula¬ 
tion  phase  and  amplitude.  Finally,  the  prognosis  for  moving  these  theoretical  and 
computational  studies  into  an  experimental  demonstration  is  commented  upon. 
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1.  PART  I  -PRINCIPLES  OF  PHOTON  MIGRATION 

1.1  Time  and  Frequency  Domain  Photon  Migration 

Photon  migration  techniques  consist  of  measuring  the  time-dependent  light  prop¬ 
agation  of  multiply  scattered  light  as  it  is  transmitted  through  tissue.  Two  techniques 
are  used  to  monitor  the  time-dependence  of  photon  migration  in  a  scattering  medium: 
(i)  time  domain  and  (ii)  frequency  domain  techniques.  In  the  time  domain  technique, 
a  picosecond  impulse  of  light  is  launched  into  a  scattering  medium  and  the  intensity  of 
the  detected  light  is  recorded  as  a  function  of  picosecond  to  nanosecond  time.  As  time 
progresses,  fewer  photons  are  detected  with  longer  “times-of-flight.”  Figure  1.2  shows 
the  broadened  pulse  that  is  re-emitted  and  measured  representing  the  distribution  of 
photon  “times-of-flight”  traveled  by  the  detected  photons.  In  the  frequency  domain 
(Figure  1.3),  incident  light  whose  intensity  is  sinusoidally  modulated  is  continuously 
launched  into  a  scattering  sample.  As  the  “photon  density  wave”  of  light  propagates, 
it  experiences  a  phase  lag  and  an  amplitude  reduction  relative  to  the  incident  light. 
The  phase-shift  and  amplitude  modulation  are  related  to  the  optical  properties  of  the 
medium.  Both  the  modulation  phase  and  amplitude  are  measured  as  a  function  of 
multiple  frequencies.  The  modulation  ratio  is  defined  as  the  modulation  amplitude 
divided  by  the  average  intensity. 

The  intensity  of  the  re-emitted  light  as  a  function  of  time  (time  domain),  or 
the  phase  and  amplitude  of  the  modulated,  re-emitted  light  as  a  function  of  the 
modulation  frequency  (frequency  domain)  contains  identical  information  regarding 
the  transport  of  light  through  the  scattering  medium.  Time  and  frequency  domain 
measurements  are  mathematically  related  to  each  other  by  the  Fourier  transform. 


6 


Photon  migration  imaging  consists  of  (i)  the  forward  problem  of  measuring  the 
time-dependent  characteristics  of  light  propagation  and  (ii)  the  inverse  problem  of 
reconstructing  an  image  of  normal  and  diseased  tissues  based  on  their  contrasting 
optical  properties.  In  the  following,  the  properties  which  impact  the  propagation  of 
excitation  light  are  reviewed. 

1.2  Transport  Property:  Absorption  Coefficient 

The  absorption  coefficient,  /j,a,  is  the  inverse  of  the  mean  distance  a  photon  will 
travel  before  being  absorbed.  In  tissue,  oxy-  and  deoxy-hemoglobin  are  the  predomi¬ 
nant  endogenous  absorbers  in  the  therapeutic  window  (Figure  1.1).  In  this  wavelength 
range,  most  tissues  have  an  absorption  coefficient  in  the  range  of  10-2  to  10-1  cm-1  [8]. 

In  the  absence  of  scattering,  the  absorption  coefficient  can  be  found  from  the 
attenuation  of  light  directly  using  the  Beer-Lambert  Law, 

Ha  =  “In  £  =  ex[C],  (1.1) 

where  ia  is  the  absorption  coefficient  of  the  medium  (cm-1);  L  is  the  optical  path 
length  light  travels  (cm);  I  is  the  intensity  of  the  detected  light  (number  of  photons/ 
(cm2-s));  I0  is  the  intensity  of  the  incident  light  (number  photons/(cm2-s));  t\  is  the 
extinction  coefficient  at  wavelength  A  (cm-1  mM'1);  and  [C]  is  the  concentration  of 
the  light  absorbing  species  (mM).  Note  that  the  number  of  photons  per  second  can 
be  converted  into  units  of  power  (Watts)  by  multiplying  by  the  photon’s  energy. 

Equation  1.1  describes  an  absorption  coefficient  that  is  inversely  related  to  the 
optical  path  length,  L.  In  the  presence  of  a  scattering  medium  such  as  tissue,  there 
is  no  unique  path  length,  but  rather  a  distribution  of  path  lengths  (Figure  1.4). 
Therefore,  determination  of  the  absorption  coefficient  in  a  scattering  medium  is  not 
possible  from  the  Beer-Lambert  relationship. 
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1.3  Transport  Property:  Scattering  Coefficient 

The  scattering  coefficient,  gs,  describes  the  distance  photons  travel  between  scat¬ 
tering  events.  Scattering  in  tissue  occurs  due  to  the  index  of  refraction  mismatch 
between  fluid  and  cellular  organelles.  Mitochondria  are  postulated  to  be  the  predom¬ 
inant  scatterer  in  tissues  [9,  10]. 

When  multiple  scattering  occurs,  light  no  longer  follows  a  straight  path.  Instead, 
photons  travel  a  distribution  of  optical  path  lengths,  rendering  the  Beer-Lambert  law 
invalid.  As  shown  in  Figure  1.4,  the  total  optical  path  length,  L,  is  now  defined  as 
the  sum  of  all  individual  path,  s,  between  scatters  and  the  scattering  coefficient  is 
defined  as  the  inverse  of  the  mean  scattering  length,  £*, 

Vs  =  ; £ ■  (1-2) 

Since  “photon  migration”  describes  a  random  walk  of  photons  between  scatter¬ 
ed,  it  is  often  more  convenient  to  define  the  scattering  parameter  as  the  isotropic 
scattering  coefficient  [//'  =  gs(l  —  g))  where  the  anisotropy  coefficient,  g,  is  defined 
below.  In  the  near  infrared  wavelength  range,  the  values  of  g's  for  tissues  typically 
range  from  10  to  50  cm-1  [8]. 

1.4  Transport  Property:  Anisotropy  Coefficient 


The  anisotropy  coefficient,  g ,  is  defined  as  the  mean  cosine  of  the  scattering  angle 
at  which  a  photon  deflects  from  a  scatterer.  Flock  et  al.  [11]  have  shown  that  the 
Henyey-Greenstein  phase  function  describes  the  angular  scatter  of  light  in  tissue.  This 
phase  function  describes  the  relationship  between  g  and  the  collimated  transmittance, 
Tc,  at  a  scattering  angle  6  where 


W)  = 


(l-<72) 


(1.3) 


(1  +  g2  —  2g  cos  0)3/2 

Collimated  transmittance  is  a  measure  of  the  light  intensity  which  is  transmitted 
through  the  sample  as  a  function  of  scattering  angle,  6.  Using  a  least  squares  analysis 
of  several  measurements  of  Tc  versus  0,  the  anisotropy  coefficient  can  be  obtained 
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from  Equation  1.3.  The  anisotropy  coefficient  varies  between  isotropic  (g  =  0)  and 
forward  (g  =  1)  scattering.  In  all  tissues,  light  is  known  to  be  forward  scattering  and 
the  anisotropy  coefficient  is  reported  to  be  around  0.9  [12].  The  Henyey-Greenstein 
phase  function  is  widely  used  in  tissue  optics  studies  not  only  because  it  provides 
a  good  representation  of  experimental  tissue  measurements  but  also  because  it  is 
mathematically  simple. 

1.5  Governing  Equations  Describing  the  Transport  of  Light 

The  scattering,  absorption,  and  anisotropy  coefficients  are  parameters  that  govern 
the  propagation  of  light  in  any  medium.  Analysis  of  light  transport  is  analogous  to 
that  of  the  motion  of  neutrons  in  a  reactor  core  where  neutrons  diffuse  from  regions  of 
high  to  low  density  [13].  When  photons  rather  than  neutrons  are  employed,  the  one- 
speed  radiative  transport  equation  describes  light  propagation  in  a  scattering  medium. 
The  radiative  transport  equation  is  written  in  terms  of  the  integro-differential  equa¬ 
tion  governing  the  local  number  of  photons,  n(r,  E,  Cl,t),  at  position  r  and  time  t, 
and  traveling  with  energy  E  and  angular  direction  Cl  [14], 

— n-~ -  =  —cCl  •  Vn  —  cgtn(r,  E,  Cl,  t)  +  s(r,  E,  Cl,  t )  + 

r  A  roo  A  A  A 

/  dCl'  /  dE'c'fi's(E'  -+  E,  Cl1  ->  fl)n(f,  E1,  Cl',  t) .  (1.4) 

The  term  on  the  left  hand  side  (LHS)  represents  the  local  accumulation  of  photons 
while  the  first  term  on  the  right  hand  side  (RHS)  represents  the  flux  of  photons  into 
the  system.  The  second  term  on  the  RHS  represents  the  loss  of  photons  traveling 
with  energy  E  in  direction  Cl  due  to  absorption  and  scattering  events.  The  total  cross 
section  for  these  extinction  processes,  /j,t,  is  given  by  the  sum  of  the  absorption  and 
scattering  cross  sections,  gt  =  fia+  g's  (cm-1).  The  third  term  on  the  RHS  represents 
the  source,  and  along  with  the  second  term  accounts  for  the  rate  of  addition  of 
photons  traveling  direction  Cl  with  energy  E  which  have  been  scattered  from  direction 
Cl'  and  energy  E' .  To  account  for  scattering  from  photons  traveling  all  directions 
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and  energies,  this  last  term  on  the  RHS  is  integrated  over  all  solid  angles,  f4w  = 
fo*  fo  sin0d0d<fr,  and  energies.  The  symbol  c  represents  the  speed  of  photons  (cm/s). 

1.5.1  One-Speed  Approximation  to  the  Transport  Equation 

When  all  the  neutrons  or  photons  travel  at  the  same  speed,  the  neutron  transport 
equation  can  be  reduced  to  the  one-speed  approximation  by  eliminating  the  depen¬ 
dence  on  energy,  E.  By  writing  the  local  concentration  of  photons  in  terms  of  an 
angular  flux,  (number  of  photons/ (cm2 •steradians-s)  where  <p(r,Cl,t)  —  cn(r,Cl,t), 
the  one-speed  equation  is  [13] 

— 7-  =  —Cl  •  V<£>  -  fit(r)<p(r,tl,t)  +  s(r,Cl,t)  +  /  dCl fJ,s(Cl -4  Cl)ip(r,Cl ,t)  .  (1.5) 

C  Ol  J 4?r 

The  term  0  •  Vy>  represents  the  negative  rate  of  addition  of  photons  traveling  an 
angular  direction  Cl.  The  term  0,  t)  accounts  for  the  loss  of  photons  traveling 

in  direction  Cl  due  to  scatter  and  absorption.  The  integral  term  represents  photons 
scattering  from  direction  Cl'  into  the  direction  Cl.  The  term  s(r,  Cl,  t)  represent  a 
source  of  photons  (number  of  photons/(cm3-steradians-s))  traveling  in  a  direction  Cl 
at  a  position  r  and  time  t  [13]. 

1.5.2  Diffusion  Equation  Describing  Light  Propagation  in  Tissues 

Following  the  derivation  of  Duderstadt  and  Hamilton  [13],  the  one-speed  equation 
can  be  reduced  to  the  diffusion  equation  describing  light  propagation  in  tissues 

~  =  V  •  D(f)V$  -  liaWQftt)  +  S(f,t ) .  (1.6) 

Equation  1.6  represents  the  time  domain  diffusion  equation  which  can  analogously 
be  written  in  the  frequency  domain  as  [14,  15] 

^Vc(f»  -  V  •  DV^c(r»  +  na$AC  =  M5(rs  -  r),  (1.7) 

where  $AC  is  the  AC  component  of  the  fluence  at  position  r  and  frequency  u>;  M  is 
the  modulation  of  the  point  source;  and  r*s  denotes  the  source  position. 
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1.5.3  Source  Location  and  Boundary  Conditions 


A  collimated  light  source  incident  on  a  scattering  medium  is  slowly  transformed 
into  a  diffuse  source  as  it  encounters  scatterers.  Patterson  et  al.  [16]  have  analytically 
shown  that  in  order  to  accurately  represent  an  incident  laser  light  source  at  the 
surface,  the  simulated  source  location  has  to  be  set  to 


(1.8) 


Therefore,  in  order  to  accurately  simulate  a  laser  source  located  on  the  surface,  the 


placement  of  the  simulated  source  needs  to  be  one  isotropic  scattering  length  inside 
the  physical  boundary.  A  sinusoidally  modulated  light  source  can  be  represented  by 
a  complex  number  [17] 

S  =  ,So[cos(0)  +  i  sin(#)] ,  (1.9) 


where  6  and  S0  are  the  phase  and  the  amplitude  of  the  source,  respectively. 

The  three  most  common  boundary  conditions  used  to  solve  the  diffusion  equation 
at  the  air-scatterer  interface  of  a  random  medium  are  the  (i)  zero  fluence  bound¬ 
ary  condition,  (ii)  the  partial  current  boundary  condition  and  (iii)  the  extrapolated 
boundary  condition.  The  zero  fluence  boundary  condition  defines  $  =  0  on  and  out¬ 
side  the  boundary.  Although  the  zero  fluence  condition  is  physically  inaccurate,  it  is 
mathematically  simple  and  gives  solutions  to  the  diffusion  equation  which  agree  well 
with  experimental  data  for  biological  systems  [8]. 


1.6  Non-Radiative  Mechanisms  for  Exogenous  Contrast 


As  discussed  above,  the  diffusion  equation  (Equations  1.6  and  1.7)  describes  the 
propagation  of  light  in  a  highly  scattering  medium.  This  equation  depends  on  three 
parameters:  (i)  the  speed  of  light,  c,  (ii)  the  isotropic  scattering  coefficient,  fi's ,  and 
(iii)  the  absorption  coefficient,  pa.  Since  all  these  parameters  govern  the  migration  of 
photons  in  tissue,  an  exogenous  compound  which  alters  the  characteristics  of  any  of 
these  parameters  can  be  used  to  induce  contrast  for  imaging  and  provide  information 
for  diagnostic  information. 
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1.6.1  Contrast  Owing  to  Absorption 

The  propagation  of  a  photon  density  wave  through  a  medium  containing  a  region 
of  light  absorbing  dye  is  different  from  the  propagation  of  a  wave  in  a  medium  with 
uniform  optical  properties  as  illustrated  in  Figure  1.5.  The  “homogeneous”  propaga¬ 
tion  is  similar  to  a  wave  formed  by  a  pebble  being  dropped  in  a  shallow  pond.  The 
propagating  wave  dampens  as  it  travels  across  the  surface  but  the  wave  maintains 
coherence.  When  the  wave  encounters  an  optical  heterogeneity  (or  diseased  tissue 
with  different  optical  properties),  it  will  reflect,  diffract,  and  be  absorbed  [18].  The 
“heterogeneous”  situation  is  similar  to  a  pebble  being  dropped  in  the  shallow  pond 
with  a  solid  obstruction  in  the  wave’s  path.  The  solid  obstruction  will  absorb  and 
partially  reflect  a  re-scattered  wave.  At  any  position,  the  diffusing  wave  will  be  the 
sum  of  the  incident  wave  and  the  scattered  wave. 

Since  these  differences  in  light  propagation  enable  contrast,  photon  migration 
measurements  are  conducted  relative  to  an  absence  (homogeneous)  condition.  The 
phase  contrast,  A0,  and  modulation  contrast,  AM,  are  defined  as 


A  6  =  6. 


presence  @ absence 


jyj  M presence 

M absence 


(1.10) 

(1.11) 


where  6presence  and  Mpresence  are  respectively,  the  measured  modulation  phase  and 
modulation  ratio  in  the  presence  of  a  heterogeneity.  Likewise,  6absence  and  Mabsence 
are  respectively,  the  measured  modulation  phase  and  modulation  ratio  in  the  absence 
of  a  heterogeneity.  If  the  optical  contrast  relative  to  an  absence  condition  is  significant 
enough  to  cause  A6  ^  0  and  AM  ^  1.0,  then  it  may  be  possible  to  reconstruct  an 
interior  optical  property  map  from  exterior  measurements  of  modulation  phase  and 
amplitude  as  described  in  Part  III. 

Since  phase  and  amplitude  contrast  depend  on  differences  caused  by  an  absorbing 
heterogeneity,  the  physics  which  change  the  detected  wave  need  to  be  examined. 
Figure  1.6  (a)  describes  frequency  domain  photon  migration  imaging  which  consists 
of  launching  a  sinusoidally  modulated  light  source  at  position  rs  and  measuring  the 
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detected  signal  at  position  fd .  The  propagation  of  the  wave  through  the  medium 
causes  a  decrease  in  modulation  and  an  increase  in  phase  at  position  r  relative  to  the 
incident  wave.  The  magnitude  of  a  scattered  wave  depends  on  the  heterogeneity’s 
distance  from  the  source,  its  dimensions,  and  the  absorption  difference  between  it  and 
the  surroundings.  The  magnitude  of  a  scattered  wave  will  be  greatest  closest  to  the 
source.  The  greatest  magnitude  possible  will  occur  when  there  is  no  loss  mechanisms 
(no  absorption)  in  the  surroundings  and  the  heterogeneity  is  a  perfect  absorber  (a 
black  body,  or  fj,a  is  infinity).  Nonetheless,  even  for  a  perfect  absorber,  the  incident 
propagated  wave  will  constitute  the  predominant  signal  measured  at  the  detector. 
Consequently,  the  scattered  wave  from  an  absorber  will  be  small  in  comparison  to 
the  incident  propagated  wave  and  the  contrast  offered  by  absorption  should  be  low. 
In  other  words,  absorption  contrast  consists  of  the  ability  to  detect  a  small  scattered 
wave  within  a  very  large  signal. 

1.6.2  Contrast  Owing  to  Dye  with  Fluorescent  Characteristics 

Another  means  to  induce  contrast  is  through  fluorescence.  Imagine  an  optical  het¬ 
erogeneity  which  contains  a  contrast  agent  that  not  only  absorbs,  but  also  fluoresces 
light  (Figure  1.6  (b)).  In  essence,  if  an  interference  filter  is  placed  at  the  detector  so 
that  only  fluorescence  light  is  seen,  the  heterogeneity  would  act  like  a  beacon.  Since 
there  is  no  small  signal  to  identify  from  a  large  background  signal,  the  contrast  should 
be  significantly  greater  than  that  offered  by  absorption.  Usually  perfect  uptake  is  not 
possible  and  one  could  imagine  that  while  there  is  10  fold  or  at  most  20  fold  more  dye 
in  a  tissue  volume  of  interest  than  in  the  surrounding  tissues,  the  problem  becomes 
one  of  picking  the  brightest  beacon  from  a  number  of  beacons  uniformly  distributed  in 
the  tissue  (Figure  1.6  (c)).  Under  these  conditions,  the  fluorescent  properties  (such  as 
lifetime  and  quantum  efficiency)  of  the  compound  may  be  used  to  impart  additional 
contrast  than  that  offered  simply  by  a  fluorophore  concentration  difference  between 
the  heterogeneity  and  its  surroundings. 
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1.6. 2.1  Physics  of  Photoluminence 

Fluorescence  is  the  emission  of  a  photon  which  results  when  a  molecule  in  an 
excited  state  relaxes  to  its  original  ground  state  (Figure  1.7).  The  excited  molecules 
can  either  relax  non-radiatively  to  the  ground  state  (without  emitting  fluorescence) 
or  stay  in  the  excited  state  for  a  period  of  time  before  returning  to  the  ground  state 
and  emitting  a  fluorescent  photon.  Often  when  pi  orbitals  are  close,  an  electron  in  the 
electronically  activated  level  will  experience  a  change  in  its  spin  state.  Since  relax¬ 
ation  to  the  ground  state  populated  with  same  spin  state  is  forbidden,  the  activated 
molecule  will  remain  in  the  excited  state  for  a  longer  period  of  time.  In  this  case  the 
emission  is  termed  phosphorescence.  Due  to  the  loss  of  energy  associated  when  the 
fluorescent  or  phosphorescent  molecules  return  to  the  ground  state,  the  emission  light 
occurs  at  a  longer  wavelength  than  the  excitation  light  [19]. 

The  fluorescent  properties  of  a  compound  are  characterized  by  the  quantum  effi¬ 
ciency,  4>i  and  the  lifetime,  r,  of  the  molecule.  The  quantum  efficiency  is  the  ratio 
of  the  number  of  fluorescent  photons  emitted  to  the  number  of  excitation  photons 
absorbed  by  the  fluorophore.  The  lifetime  of  the  molecule  is  the  average  time  that 
the  molecule  stays  in  the  activated  state  and  is  measured  as  the  average  time  be¬ 
tween  absorption  of  an  excitation  photon  and  emission  of  a  fluorescent  photon.  Both 
lifetime  and  quantum  efficiency  are  sensitive  to  the  local  bio-chemical  environment. 
For  example,  porphyrin  fluorescence  is  quenched  by  oxygen  enabling  quantitative  de¬ 
tection  of  oxygen  from  the  shortening  of  its  lifetime  and  reduction  of  its  quantum 
efficiency  [19]. 

1.6. 2. 2  Measurement  of  Photolumination  Properties  in  Dilute,  Non-scatter¬ 
ing  Suspensions 

The  lifetime,  r,  and  quantum  efficiency,  <f>,  of  a  molecule  can  be  determined  using 
both  a  time  and  frequency  domain  analysis.  In  the  time  domain,  a  pulse  of  light 
excites  a  dilute  non-scattering  solution  of  fluorophores  and  the  decay  of  the  gener¬ 
ated  fluorescent  intensity  is  recorded  as  a  function  of  time  after  the  initial  impulse. 
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The  amount  of  fluorescent  light,  /m,  is  defined  by  the  absorption  coefficient  of  the 
fluorophores,  fiaxf,  the  quantum  efficiency,  and  the  lifetime,  r,  of  the  molecule  [19], 

Im(t)  =  ttlzLe~i/T  _  I0e~l!\  (1.12) 

T 

where  I0  is  the  maximum  fluorescent  intensity.  The  lifetime  of  a  fluorophore  is  equal 
to  the  time  required  for  the  intensity  to  decay  to  1/e  of  its  initial  value. 

Frequency  domain  measurements  of  lifetime,  r,  and  quantum  efficiency,  cf>,  are 
made  by  exciting  a  dilute,  non-scattering  solution  of  fluorophore  with  a  sinusoidally 
modulated  source  and  measuring  either  the  emitted  modulation  phase  or  modulation 
ratio  of  the  fluorescent  signal  relative  to  the  incident  excitation  light.  The  modulation 
phase  of  the  fluorescent  light  is  independent  of  the  fluorophore  concentration  and  the 
quantum  efficiency, 

;  0m(u>)  =  tan-^wr),  (1.13) 


where  u  is  the  modulation  frequency  of  the  light  source.  The  modulation  ratio, 
however,  is  dependent  on  quantum  efficiency,  </>,  and  [iaxf , 


Mm(u>) 


l/1  + 


(1.14) 


and  can  therefore  be  used  to  determine  both  lifetime,  r,  and  quantum  efficiency,  <f>,  of 
a  sample.  If  the  modulation  ratio  is  referenced  to  the  modulation  ratio  obtained  under 
continuous  wave  or  constant  illumination  (u>  =  0),  the  resulting  ratio  is  a  function  of 
modulation  frequency  and  lifetime,  r, 


Mm(u)  _  1 

Mm( 0)  “  +  (wt)2  ' 


(1.15) 


Equations  1.13  -  1.15  are  obtained  from  the  Fourier  transform  of  Equation  1.12. 

The  measurement  of  the  emission  properties  of  r  and  <f)  are  determined  in  a  dilute, 
non-scattering  medium.  Since  the  fluorescent  lifetime,  r,  is  on  the  same  order  as 
photon  “times-of-flight,”  the  equations  above  are  not  valid  in  a  scattering  medium  due 
to  the  additional  time  delay  and  amplitude  decrease  associated  with  photon  migration. 
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In  order  to  accurately  determine  the  emission  properties  in  a  scattering  medium, 
photon  migration  times  need  to  be  deconvolved  from  time  dependent  measurements. 

As  an  alternative,  long  lived  (phosphorescent)  probes  have  been  suggested  for 
photon  migration  measurements  of  the  re-emission  properties  in  a  highly  scattering 
medium  since  the  phosphorescent  lifetimes  are  longer  than  photon  “times-of-flight.” 
However,  Hutchinson  et  al.  [20,  21]  found  that  the  signal  originating  from  a  medium 
containing  phosphorescent  agents  is  confined  to  surface  or  sub-surface  regions.  Since 
phosphorescent  agents  cannot  be  used  to  determine  the  lifetime  deep  within  scattering 
media,  fluorescent  contrast  agents  are  suggested. 

1.6. 2. 3  Collisional  Quenching  Mechanisms 

Collisional  quenching  is  the  mechanism  in  which  the  fluorescent  intensity  decreases 
due  to  collisional  encounters  between  a  fluorophore  and  a  quencher.  For  this  process 
to  happen,  the  quencher  must  diffuse  to  the  fluorophore  while  the  fluorophore  is 
in  the  excited  state.  Some  examples  of  collisional  quenchers  are  oxygen,  chloride, 
chlorinated  hydrocarbons,  xenon,  hydrogen  peroxide,  bromate,  and  iodide  [19].  The 
lifetime  of  a  fluorophore  can  be  related  to  the  concentration  of  a  quencher  using  the 
Stern- Volmer  equation  [19], 

-  =  l  +  KSv[Q],  (1-16) 

T 

where  r0  is  the  lifetime  of  the  probe  in  the  absence  of  a  quencher,  r  is  the  lifetime  of  the 
probe  in  the  presence  of  a  quencher,  I<sv  is  the  Stern- Volmer  constant  and  [Q]  is  the 
quencher  concentration.  As  the  concentration  of  the  quencher  increases,  the  excited 
state  of  the  probe  is  quenched  causing  a  reduction  in  the  lifetime.  Since  the  (bio-) 
chemical  environment  (amount  of  quencher  present)  impacts  fluorescence  lifetime,  r, 
fluorescent  lifetime  imaging  may  be  used  for  metabolite  sensing.  It  is  important  to 
note  that  Equation  1.16  does  not  depend  on  the  concentration  of  fluorophore. 
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A  class  of  lifetime  sensitive  fluorophores,  called  multi-plex  dyes  [22]  involve  “tun¬ 
ing”  cyanine  dyes  to  exhibit  different  spectra  lifetimes  depending  on  local  environ¬ 
ments  and  binding  characteristics.  This  ability  to  change  the  emission  characteris¬ 
tics  between  normal  and  diseased  tissue  using  dyes  whose  characteristics  change  by 
varying  pH  [1,  2]  and  O2  [3]  may  not  only  enhance  detection  but  also  may  contain 
diagnostic  information  regarding  the  disease. 


1.6.3  Diffusion  Equations  Describing  Excitation  and  Emission  Photon 
Migration 

To  use  fluorescent  contrast  agents  in  the  forward  and  the  inverse  imaging  prob¬ 
lems,  fluorescent  generation  and  propagation  must  be  accurately  modeled.  The  light 
propagation  due  to  fluorescence  in  a  scattering  medium  can  be  described  by  the  cou¬ 
pled  diffusion  equations  for  light  propagation  [20,  23,  24], 

r  ioj 


-V  •  Z)*(r)V$*(r»  + 


—  +  Vajr)  +  Vaxf{r) 
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-V  •  Dm(f)V$m(r, w)  +  \—+fiam(r) }  $»)  =  fra- 

-Cm  J  * 
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(1.17) 

+  (L1S) 

In  the  above  equations,  <hr  and  are  the  AC  components  of  fluence  for  excitation 

(superscript  x )  and  emission  (superscript  m)  light  (photons/cm2);  /u,axi  is  the  absorp¬ 
tion  due  to  chromophores  (cm-1);  fiaxf  is  the  absorption  due  to  fluorophores  (cm-1); 
Ham  represent  the  absorption  of  the  emission  light  due  to  chromophores  (cm-1);  fiSx 
and  /jSm  are  the  scattering  coefficients  of  excitation  and  emission  light  (cm-1),  re¬ 
spectively;  4>  and  r  are  the  quantum  efficiency  and  lifetime  (ns)  of  the  fluorophore, 
respectively;  and  Dx  and  Dm  are  the  optical  diffusion  coefficients  for  the  excitation 
and  emission  light  (cm)  where 

1 


Dx  3[//0l  +  (1  -  g)pu] 


and 


Dn 


(1.19) 


(1.20) 


Hp  +  (1  -9)1*  «Sm  ) 

Equations  1.17,  and  1.18  are  coupled  by  the  excitation  fluence,  $*,  and  the  absorption 
due  to  fluorescence,  HaxJ-  The  solutions  for  these  equations,  and  $m,  are  complex. 
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It  is  important  to  note  that  fluorescence  is  assumed  linear  with  the  excitation  fluence 
(e.g.  photo-bleaching  invalidates  the  model).  The  modulation  phase,  6 ,  and  modula¬ 
tion  amplitude,  Iac ,  of  the  excitation  and  emission  light  are  represented  by  the  real 
and  imaginary  components  of  the  emission  and  excitation  fluence, 


ex’m  =  tan-1 


'Im§x'm' 
.  Re$x’m . 


(1.21) 


(1.22) 


and  the  modulation  ratio  is  obtained  by  dividing  Iac  by  the  average  intensity  (or  Idc), 


I  DC 


(1.23) 


The  mathematical  framework  for  predicting  modulation  phase  and  modulation 
amplitude  at  emission  wavelengths  assumes  an  angle  exponential  decay  kinetics,  but 
can  be  easily  modified  for  other  kinetics. 


1.6.4  Potential  Fluorescent  Contrast  Agents  for  Photon  Migration  Imag¬ 
ing 

An  example  of  a  class  of  potential  optical  contrast  agents  are  porphyrin  derivatives. 
Porphyrin  molecules  have  characteristic  ring  structures  derived  from  four  pyrrole  rings 
joined  together  by  four  methene  bridges  as  shown  in  Figure  1.8.  Porphyrin  derivatives 
are  typically  employed  in  photodynamic  therapy  (PDT).  PDT  is  a  method  in  which 
a  hematoporphyrin  derivative  drug  is  administered  for  photo-therapeutic  destruction 
of  cancer  [25].  In  this  therapy,  a  photosensitizing  drug,  which  has  a  high  affinity 
for  tumor  tissue,  is  systematically  administered  to  the  patient.  The  tumor  area 
is  then  irradiated  with  red  light  which  photoactivates  the  dye  to  a  triplet  state  to 
produce  cytotoxic  by-products  causing  irreversible  cellular  damage.  Many  problems 
limit  the  efficacy  of  photodynamic  therapy  such  as  damage  to  skin  exposed  to  sunlight. 
Nonetheless,  while  PDT  agents  depend  upon  triplet  states  for  their  therapeutic  action, 
the  singlet  state  may  enable  contrast  for  diagnostic  imaging. 
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Indocyanine  Green  (ICG)  (Figure  1.9)  is  another  optically  active  contrast  agent 
that  is  approved  by  the  Food  and  Drug  Administration  (FDA)  for  diagnostic  studies. 
It  is  currently  used  to  diagnose  retinal  and  chorodial  diseases  by  enhancing  the  images 
of  the  retinal  vasculature  of  the  eye  [26].  In  burn  diagnostics,  ICG  is  administered  to 
the  blood  stream  after  which  the  burn  area  is  illuminated  to  monitor  the  amount  of 
blood  flow  to  the  area.  Because  blood  flow  is  proportional  to  the  depth  of  the  burn, 
ICG  is  used  to  measure  the  severity  of  the  wound  [27].  ICG  is  also  used  to  determine 
hepatic  function  by  monitoring  the  amount  of  dye  that  clears  the  circulatory  system 
of  the  liver  as  function  of  time.  Sevick  et  al.  [4]  measured  the  extinction  coefficient, 
ex,  of  ICG  to  be  1,300  (M  cm)-1  at  780  nm  and  22,000  (M  cm)-1  at  830  nm.  They 
also  report  the  quantum  efficiency,  4>,  and  lifetime,  r,  of  ICG  to  be  0.016  and  0.56  ns, 
respectively.  Although  their  measurements  were  conducted  with  ICG  diluted  in  water, 
the  properties  are  assumed  to  be  reflective  of  those  in  0.5%  Intralipid. 
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2.  PART  II-CONTRAST  FOR  PHOTON  MIGRATION 
IMAGING:  FLUORESCENCE  LIFETIME 

2.1  Single  Pixel  Measurements  for  Photon  Migration  Imaging 

To  experimentally  assess  the  feasibility  of  external  contrast  agents,  the  impact  of 
absorbing  and  fluorescing  heterogeneities  on  light  propagation  in  a  tissue  simulating 
phantom  was  examined.  Also,  the  influence  of  lifetime  on  photon  migration  mea¬ 
surements  was  investigated.  Single  pixel  measurements  of  modulation  phase,  0,  and 
modulation  ratio,  M,  were  measured  as  a  function  of  modulation  frequency,  w,  in  the 
presence  and  absence  of:  (i)  light-absorbing,  and  (ii)  fluorescing  objects. 

2.2  Phantom  Apparatus 

The  tissue  phantom,  illustrated  in  Figure  2.1,  consisted  of  a  Plexiglas  cylindrical 
tank,  20  cm  in  diameter  and  30.5  cm  in  height,  filled  with  8  liters  of  0.5%  by  volume 
Intralipid  solution  (Kabi  Pharmacia,  Clayton,  NC).  This  concentration  of  Intralipid 
mimics  the  scattering  properties  of  tissue  (//'  «  10  cm-1)  [28].  Two  small  holes,  an 
arc  distance  of  2.0  cm  apart,  were  drilled  half  way  up  the  side  of  the  tank  in  order  to 
couple  the  1000  micron  fiber  optic  (HCP-M1000T-08,  Spectron  Specialty  Optics  Co., 
Avon,  CT)  source  and  detector  directly  to  the  phantom.  The  fiber  optics  were  placed 
approximately  1  mm  inside  the  medium  to  avoid  boundary  effects  caused  by  the  wall 
of  the  tank.  On  the  outside  of  the  tank,  silicon  gel  (Dow  Corning  Co.,  high  vacuum 
grease,  Midland,  MI)  was  placed  around  the  fibers,  to  prevent  leakage  of  Intralipid. 

The  perfect  absorbing  heterogeneity  was  a  6  mm  diameter  glass  rod  painted 
black.  The  fluorescent  body  consisted  of  a  9  mm  inner  diameter  cylindrical  glass 
container  filled  with  micromolar  solutions  of  either  3,3’-Diethyl-thiatricarbocyanine 
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iodide  (DTTCI)  (ACROS  Organics,  Fairlawn,  NJ),  a  common  laser  dye,  or  ICG 
(Sigma  Chemical  Company,  St.  Louis,  MO)  diluted  in  0.5%  Intralipid  solution.  Di¬ 
lution  in  Intralipid  maintained  a  scattering  coefficient  inside  the  heterogeneity  to  be 
identical  to  the  surrounding  medium.  Measurements  of  6  and  M  were  conducted  as 
the  heterogeneity  was  moved  various  distances  away  from  the  wall  midpoint  between 
the  source  and  detector.  The  object  position  was  measured  from  the  leading  edge  of 
the  heterogeneity  where  a  position  of  zero  denotes  contact  between  the  heterogeneity 
and  the  wall  of  the  tank. 

2.3  Frequency  Domain  Instrumentation  and  Set-Up 

The  light  propagation  at  both  the  excitation  and  emission  wavelengths  was  mea¬ 
sured  using  frequency  domain  measurements.  The  measurements  were  conducted  at 
modulation  frequencies  between  80  and  240  MHz  using  a  pulsed  Ti:Sapphire  laser, 
coupled  to  the  detection  system  consisting  of  an  ISS  K2  phase  fluorometer  (ISS, 
Champaign,  IL). 

2.3.1  Heterodyne  Detection  Approach  with  Modulated  Laser  Source 

To  use  simple  electronics  to  detect  the  signal,  a  heterodyne  technique  was  used. 
This  detection  method  involves  modulating  the  detectors  at  a  frequency  plus  a  small 
offset  from  the  modulation  frequency  of  the  source.  For  example,  the  sinusoidal  source 
is  described  by  the  expression 

A  cos (u>t  +  0),  (2-1) 

where  A  is  the  amplitude,  u>  is  the  modulation  frequency  and  9  is  the  phase  of  the 
source.  The  gain  of  the  detector  is  modulated  by  w+Aw  so  that  the  detector  response 
is  represented  by 


cos[(u>  +  Au>)t], 


(2.2) 
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where  Atu  is  the  small  offset  frequency  of  the  detector.  Then  the  signal  obtained  by 
the  detector  is 

A  A 

Acos(u>t  +  0)  x  cos[(u;  +  Aa>)f]  =  —  cos(A ut  +  0)  +  —  cos[(2tu  +  A ui)t  +  0].  (2.3) 

A  low  pass  filter  is  used  to  filter  out  frequencies  above  Atu,  the  signal  received  by  the 
detector  is 

—■  cos(A wt  +  6).  (2.4) 

A 

It  is  important  to  note  that  the  phase  and  amplitude  information  of  the  high  frequency 
signal  is  preserved. 

2.3.2  Heterodyne  Detection  Approach  with  Pulsed  Laser  Source 

When  the  light  source  is  a  mode  locked  pulsed  laser,  the  signal  can  no  longer  be 
described  by  Equation  2.1.  Instead,  the  source  is  represented  by 

OO 

^2  An  cos  (nut  -1-  0n),  (2.5) 

71=1 

where  the  above  equation  is  the  frequency  domain  representation  of  a  mode  locked 
pulsed  laser  [29].  Again  the  detector  is  modulated  at  a  frequency  plus  a  small  offset 
as  described  by  Equation  2.2,  the  signal  obtained  by  the  detector  is 

cos[(tu  +  Acu)t]  x  [Ai  cos(u>t  +  6i)  +  A 2  cos(2cut  +  62)  +  . . .]  = 

T  ~z~ "{cos(Atuf  $1)  T  cos[(2tu  T  At u)t  T  $1]}  (2.6) 

Ao 

+  -^-{cos[(tu  +  At S)t  +  $2 ]  +  cos [(3a;  +  At o)t  +  ^2]}  T  •  •  • 
Filtering  out  frequencies  above  At u,  the  signal  received  by  the  detector  is 

— -  cos(  At vt  +  61 )  (2.7) 

which  is  completely  analogous  to  Equation  2.4. 

2.3.3  Instrumentation 

The  excitation  light  source  consists  of  a  Ti: Sapphire  laser  pumped  by  a  10  Watt 
Argon  ion  laser  (Spectra  Physics,  CA)  to  produce  a  pulse  train  of  2  ps  at  a  repetition 
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rate  of  80  MHz  (Figure  2.2).  The  output  wavelength  of  the  Ti:Sapphire  laser  is  set 
to  780  nm.  To  establish  a  phase  lock  with  the  detection  equipment,  an  80  MHz 
signal  origination  from  the  Ti:Sapphire  laser  is  sent  to  a  divide-by-8  circuit  resulting 
in  a  10  MHz  signal.  This  10  MHz  signal  is  sent  to  a  frequency  synthesizer  (Model 
2022D,  Marconi,  Allendale,  NJ)  which  controls  the  modulation  of  the  detectors.  The 
detectors  are  modulated  at  frequencies  of  80,  160  and  240  MHz.  The  light  from  the 
TkSapphire  laser  has  1.3  Watts  of  average  power. 

The  pulse  train  is  focused  onto  a  beam  splitter.  The  unreflected  portion  of  the 
pulse  train  (approximately  90%  of  the  signal)  serves  as  the  light  source  and  is  focused 
onto  a  fiber  optic  coupled  to  the  side  of  the  cylindrical  tank.  The  reflected  portion 
of  the  signal  is  focused  onto  another  fiber  optic  cable  which  delivers  the  signal  to 
a  reference  photomultiplier  tube  (PMT).  This  light  serves  as  the  reference  signal. 
The  scattered  excitation  and  emission  light  detected  at  the  periphery  of  the  phantom 
is  delivered  to  a  sample  PMT  also  using  a  fiber  optic  cable.  An  830  nm  interfer¬ 
ence  filter  (10-830-4-1.00,  CVI  Laser  Co,  Albuquerque,  NM)  is  placed  in  front  of  the 
sample  PMT  and  is  used  to  make  frequency  domain  measurements  at  the  emission 
wavelengths. 

The  gain  of  the  PMTs  are  modulated  at  a  harmonic  of  the  source  frequency  but 
offset  by  100  Hz  in  order  to  produce  a  heterodyne  signal.  The  output  from  both  PMTs 
is  sent  to  an  ISS  board  for  data  acquisition.  The  modulated  phase  and  amplitude  are 
extracted  from  the  signal  at  the  100  Hz  cross  correlation  frequency  through  simple 
electronics  and  software  (ISS,  Champaign,  IL). 

2.4  Contrast  Owing  to  Perfect  Absorption  and  Perfect  Uptake  of  Fluo¬ 
rescence 

The  modulation  phase  and  modulation  ratio  were  measured  to  evaluate  the  con¬ 
trast  offered  by  (i)  a  perfect  absorber  and  (ii)  a  fluorescent  volume  in  a  lossless 
scattering  medium  (no  absorption  in  the  surrounding).  Since  the  fluorescent  object 
had  perfect  uptake  of  dye  and  since  the  absorber  was  a  black  rod,  they  represent 
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the  best  possible  contrast  owing  to  absorption  and  fluorescence.  For  this  experiment 
DTTCI  served  as  the  fluorescent  contrast  agent. 

Measurements  of  modulation  phase,  9,  and  modulation  ratio,  M,  were  conducted 
as  the  heterogeneities  were  moved  from  the  peripheral  location  towards  the  center 
of  the  phantom.  Figure  2.3  shows  the  phase  contrast,  A 6  ( 9presence  —  9 absence),  due 
to  the  presence  of  both  an  absorbing  and  a  fluorescing  heterogeneity  at  80,  160,  and 
240  MHz  modulation  frequencies.  The  largest  phase  contrast  is  measured  at  240 
MHz  and  the  signal  is  only  altered  by  approximately  20°  for  absorption,  however  for 
fluorescence,  an  alteration  of  approximately  60°  is  observed.  Therefore,  these  results 
show  that  the  best  contrast  available  from  absorption  (i.e.,  a  perfect  absorber)  causes  a 
smaller  measurable  disturbance  in  photon  migration  than  a  micromolar  concentration 
of  fluorescent  dye.  Also,  the  fluorescent  signal  is  detected  approximately  1  cm  deeper 
than  the  absorbance  signal. 

Measurements  of  modulation  ratio,  M,  for  the  perfect  absorber  and  the  fluorescent 
volume  were  not  performed.  However,  a  similar  experiment  conducted  by  Lopez  [30] 
using  a  9  mm  diameter  perfect  absorbing  fluorescent  object  shows  that  the  modulation 
contrast  ( Mpresence / Mabsence )  of  the  fluorescent  object  at  240  MHz  is  approximately 
0.1  while  the  absorber  only  offers  0.6  contrast.  These  results  are  consistent  with  what 
is  expected. 

In  a  realistic  situation,  the  administration  of  contrast  agents  results  in  the  imper¬ 
fect  uptake  of  dye  into  tissue  regions  of  interest.  Recently,  Lopez  et  al.  [4,  30]  have 
conducted  measurements  of  phase  and  modulation  contrast  with  imperfect  uptake 
ratios  of  1:100  and  1:10  with  a  similar  experimental  apparatus  as  described  above. 
They  found  the  resulting  phase  contrast  to  be  dramatically  smaller  than  the  situa¬ 
tion  of  perfect  uptake.  However,  the  magnitude  of  contrast  offered  by  fluorescence 
phase  contrast  and  modulation  contrast  still  exceeds  absorption  contrast,  which  is  in 
agreement  with  the  work  above.  As  discussed  in  Part  I,  the  additional  time-delay  as¬ 
sociated  with  the  lifetime  of  a  fluorescent  molecule  during  the  emission  process  causes 
a  greater  phase  and  modulation  contrast  over  that  of  absorption. 
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2.5  Contrast  Owing  to  the  Influence  of  Lifetime 

The  next  experiment  assessed  the  influence  of  lifetime  on  contrast.  In  an  analogous 
single  pixel  experiment  to  the  one  described  above,  Intralipid  solutions  with  micromo¬ 
lar  concentrations  of  DTTCI  and  ICG  were  investigated  and  compared  to  each  other 
since  DTTCI  and  ICG  have  different  lifetimes  of  1.18  and  0.56  ns,  respectively  [4]. 

Dilute  non-scattering  solutions  of  both  DTTCI  and  ICG  were  separately  analyzed 
in  a  small  glass  cuvette.  Using  the  instrumentation  discussed  above,  measurements 
of  modulation  phase  and  modulation  ratio  were  conducted  with  the  source  at  a  right 
angle  from  the  detector.  Using  a  third  dye,  IR-140  (ACROS  Organics,  Fairlawn, 
NJ),  as  a  standard,  the  lifetimes  of  DTTCI  and  ICG  were  obtained  by  evaluating  the 
equations  discussed  in  Section  1.6. 2.1. 

Both  the  modulation  phase  and  modulation  ratio  were  measured  as  a  function  of 
frequency  as  DTTCI  and  ICG  scattering  solutions  were  moved  from  the  periphery  to 
the  center  of  the  phantom.  Figure  2.4  and  Figure  2.5  illustrate  the  phase  contrast 
for  DTTCI  and  ICG,  respectively.  As  the  lifetime  of  the  dye  increases,  the  phase 
contrast  should  also  increase  and  the  modulation  ratio  should  decrease  (see  Part  I). 
Since  the  lifetimes  of  DTTCI  and  ICG  are  1.05  and  0.58  ns  respectively,  the  phase 
contrast  for  DTTCI  should  be  larger  and  the  modulation  contrast  should  be  smaller. 
This  is  indeed  the  case  since  DTTCI  shows  values  of  A 0  =  60°  and  AM  =  0.3  and 
ICG  shows  values  of  A0  =  46°  and  AM  =  0.5.  These  results  confirm  that  frequency 
domain  measurements  contain  information  of  fluorescent  lifetime. 

2.6  Discussion 

The  results  in  this  section  show  that  fluorescence  offers  improved  contrast  over 
absorption.  Since  measurements  were  conducted  on  a  perfect  absorber,  the  findings 
suggest  that  the  endogenous  contrast  due  to  absorption  from  increased  hemoglobin 
associated  with  angiogenesis  may  be  too  small  for  realistic  measurements  for  optical 
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tomography  (Figure  1.6).  The  second  study  shows  that  frequency  domain  measure¬ 
ments  contain  lifetime  information  which  can  provide  diagnostic  information  about 
the  local  bio-chemical  environment  (see  Equation  1.16).  However,  these  measure¬ 
ments  need  to  be  coupled  to  an  inversion  algorithm  in  order  to  extract  the  lifetime 
from  the  scattered  signal  (see  Part  III  for  more  details).  The  single  pixel  measure¬ 
ments  presented  here  show  great  promise  for  the  detection  of  diseased  tissues  laden 
with  contrast  agents  using  photon  migration  measurements.  But  single  pixel  mea¬ 
surements  only  provide  one  dimensional  information  and  therefore  cannot  be  used 
to  resolve  an  interior  image  of  optical  properties.  To  obtain  enough  information  to 
locate  a  heterogeneity,  measurements  at  multiple  locations  need  to  be  obtained. 
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3.  PART  III-RECONSTRUCTION  OF  FLUORESCENCE 
QUANTUM  EFFICIENCY  AND  LIFETIME 

3.1  INVERSE  IMAGING  PROBLEM 

Recently,  we  have  developed  multi-pixel  measurements  for  rapid  data  acquisition 
of  modulation  phase  and  amplitude  of  re-emitted  optical  signals  [31].  Using  a  gain- 
modulated  image  intensifier  coupled  to  a  CCD  camera,  acquisition  of  128  X  128 
measurements  of  modulation  phase  and  amplitude  are  enabled  in  approximately  10 
milliseconds  at  excitation  wavelengths  and  1  second  at  emission  wavelengths.  In  order 
to  employ  these  measurements  to  reconstruct  interior  maps  of  lifetime  and  quantum 
efficiency  from  exterior  frequency  domain  measurements  made  at  the  periphery,  an 
inverse  algorithm  needs  to  be  employed.  Currently,  two  approaches  for  solving  the 
inverse  imaging  problem  have  been  demonstrated.  Both  approaches  involve  solutions 
to  the  diffusion  equation  where  the  first  method  uses  perturbation  expansions  of  the 
photon  diffusion  equation  [32,  33,  34,  35]  and  the  second  solves  the  diffusing  equation 
directly  using  a  finite  difference  or  a  finite  element  approach  [36,  37,  38,  39,  40]. 

In  this  work,  we  employ  the  latter  approach  to  reconstruct  optical  property  maps 
of  absorption,  fluorescent  quantum  efficiency,  and  fluorescent  lifetime  of  an  optical 
heterogeneity  without  a  priori  information  of  position  or  volume.  Specifically,  the 
tissue  volume  of  interest  is  mathematically  represented  by  discretized  volume  ele¬ 
ments  in  which  the  optical  properties  of  absorption,  fluorescence  quantum  efficiency 
and  lifetime  are  constructed.  The  inversion  schemes  incorporate  an  iterative  Newton- 
Raphson  technique  to  update  values  of  the  optical  properties  at  every  volume  element 
in  order  to  minimize  the  least-squares-difference  between  the  experimental  and  cal¬ 
culated  detector  responses  at  the  periphery.  The  method  is  similar  to  that  used  by 
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Yorkey  et  al.  [41]  for  electrical  impedance  tomography,  Pogue  et  al.  [37]  for  ab¬ 
sorption  and  scattering  reconstructions,  and  Paithankar  et  al.  [40]  who  performed 
reconstructions  on  both  fluorescent  lifetime,  r,  and  fluorescent  yield  (the  product 
of  the  absorption  coefficient  due  to  fluorophores  and  quantum  efficiency,  4>/j,axf).  In 
this  section,  interested  readers  will  find  the  numerical  technique  used  in  solving  the 
forward  problem  along  with  the  approach  to  solving  the  inverse  imaging  problem. 
Finally,  applications  of  the  inversion  to  simulated  data  are  presented  in  Part  IV  illus¬ 
trating  the  feasibility  of  fluorescent  quantum  efficiency  and  lifetime  imaging  in  tissues 
and  other  scattering  media.  Part  III  may  be  omitted  by  those  readers  un-interested 
in  the  mathematics  of  solving  the  forward  and  inverse  imaging  problems. 

3.1.1  Forward  Problem 

The  coupled  diffusion  equations  for  excitation  and  emission  light  propagation 
(Equations  1.17  and  1.18)  are  complex  linear  elliptic  equations  which  can  be  analyti¬ 
cally  solved  as  boundary  value  problems  when  the  optical  properties  are  independent 
of  spatial  position.  The  complex  fluence  can  be  used  to  predict  modulation  phase 
and  amplitude  from  the  real  and  imaginary  components.  However,  when  there  is  a 
spatial  variation  in  optical  properties  (such  is  the  case  when  an  optical  heterogene¬ 
ity  is  present)  numerical  solutions  are  necessary.  In  this  work,  the  finite  difference 
method  was  used  to  solve  for  both  the  excitation  fluence,  $x(r,uj),  and  the  emission 
fluence,  4>m(r, u).  The  method  of  finite  differences  involves  discretizing  the  area  of 
interest  into  a  grid  and  a  linear  system  is  obtained  by  approximating  derivatives  in 
the  diffusion  equation  by  differences. 

In  this  study,  a  two-dimensional  simulated  tissue  phantom  as  shown  in  Figure  3.1 
was  employed.  The  phantom  is  modeled  as  a  4  x  4  cm2  square  surface  with  a  17  x 
17  grid  representing  the  area  of  potentially  differing  optical  properties.  The  source 
and  the  detector  locations  were  placed  one  pixel  in  from  the  edge  of  the  phantom  to 
overcome  the  effects  of  the  zero  fluence  boundary  condition.  A  source  is  centered  on 
each  side  of  the  phantom  for  a  total  of  4  sources.  The  grid  contains  56  detectors  per 
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source  evenly  spaced  around  the  phantom.  The  placement  of  the  source  represents  an 
isotropic  excitation  source  located  0.25  cm  or  2.5  scattering  lengths  from  the  surface. 
Since  all  the  experiments  discussed  here  are  performed  using  simulated  data,  the 
placement  of  the  source  will  not  affect  the  solution  to  the  inversion.  However,  when 
this  algorithm  is  coupled  to  experimental  measurements,  the  placement  of  the  source 
needs  to  be  correctly  represented. 

For  solutions  of  the  forward  problem,  the  known  optical  properties  of  the  back¬ 
ground  and  the  object  are  input  as  parameters  into  the  program  and  then  the  fluence 
at  each  grid  point  is  calculated.  The  fluence  values  are  used  at  the  detector  positions 
and  the  modulation  phase,  0x,m,  and  the  log  of  the  modulation  amplitude,  I^ci  are 
calculated  from  both  the  real  and  imaginary  components  of  the  excitation  or  emission 
fluence,  where 


6x'm  =  tan-1 


(3.1) 


I%c  =  log10  .  (3.2) 

To  run  a  single  simulated  experiment,  only  one  source  is  active.  For  each  forward 
solution,  values  of  9x,m  and  I^'c  are  calculated  at  all  the  detector  positions.  Once 
these  detector  values  are  obtained,  the  forward  solution  is  repeated  for  the  next 
source  location.  This  procedure  continues  until  all  four  sources  are  investigated.  A 
total  of  four  simulated  experiments  constitute  the  data  to  be  acquired  experimentally. 
Consequently,  the  simulations  will  yield  a  total  of  224  detector  measurements.  In  all 
the  simulated  experiments,  the  sources  are  modulated  at  a  frequency  of  100  MHz. 
Typical  forward  calculations  on  a  Ultra  2  Sun  workstation  took  approximately  1.3  s. 

After  detector  data  is  evaluated  for  6x,m  and  7^’™,  simulated  Gaussian  noise  with 
a  standard  deviation  of  0.1°  for  phase  and  1%  for  modulation  amplitude  is  added 
using  a  MATLAB  routine.  The  Gaussian  noise  is  added  to  account  for  the  noise 
which  would  be  encountered  in  a  real  experiment.  The  simulated  forward  results 
with  added  noise  are  then  input  to  the  inversion  algorithm. 
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3.1.2  Solution  to  the  Inverse  Problem 


The  solution  to  the  inverse  problem  requires  reconstructing  images  of  the  interior 
volume  from  exterior  measurements  of  0x,m  and  IAq.  Since  more  information  from 
exterior  measurements  enable  better  image  reconstruction,  more  detectors  and  mul¬ 
tiple  modulation  frequencies  should  improve  image  reconstructions.  The  inversion 
algorithm  is  conducted  in  two  parts  on  the  basis  of  (i)  absorption  at  the  excitation 
wavelength  due  to  fluorophores,  iiaxf ,  and  (ii)  quantum  efficiency,  <f>,  or  lifetime,  r,  at 
the  emission  wavelength.  The  two  parts  are  similar  except  the  reconstruction  based 
on  naif  uses  detector  data  at  the  excitation  wavelength  while  the  other  part  uses  data 
at  the  emission  wavelength. 

Figure  3.2  depicts  the  flow  diagram  for  these  inversion  algorithms.  To  begin 
the  inverse  calculation  based  on  Haxf,  an  initial  homogeneous  guess  for  the  optical 
property  map  is  given.  Typical  values  are  those  found  in  normal  tissues.  From 
these  values  the  forward  solution  is  calculated  to  obtain  6X  and  IXAC  measurements 
at  the  detector  positions  for  the  excitation  light.  Next,  two  Jacobian  matrices, 
J(9x,lj,axf)  and  are  constructed  where  an  entry  of  each  matrix  is  de¬ 

fined  as  jij  =  d9f /d/j,axfj  and  jtj  =  dIAC./diiaxf,j,  respectively.  Each  element  of 
the  matrix  describes  the  response  of  the  source-detector  pair  at  position  i  to  changes 
in  fiax}  at  each  pixel  j.  The  partial  derivatives  are  numerically  approximated  by 
calling  the  forward  solution  (Cud2)  a  second  time  for  a  1%  increase  of  the  original 
pixel  value  for  piaxf  and  subtracting  the  difference  (i.e.  dOx/dfiaxf  «  AOx/Afj,Uxf  and 
dIXAc/d»axf  *  AIac/ A/J,axf). 

To  update  an  iterative  Newton-Raphson  technique  involving  the  Jacobians 
was  used, 
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where  A ftax/  provides  an  increment  change  for  updating  fiaxf  and  £  is  a  regularization 
parameter  multiplied  by  the  identity  matrix,  I.  Because  the  Jacobian  matrices  are  ill 
conditioned  due  to  the  small  sensitivity  of  far  away  from  the  source  and  detector, 
the  regularization  parameter  compensates  by  making  the  matrices  more  diagonally 
dominant.  The  parameter  £  is  then  adjusted  by  a  Marquardt-Levenberg  algorithm  at 
every  iteration  [37].  In  order  to  solve  the  linear  algebraic  equations  (Equation  3.3)  for 
A jlaxf ,  a  LU  decomposition  back  substitution  method  is  used.  The  updated  optical 
property  map  is  then  found  by  adding  A jxaxf  to  the  values  from  the  previous  iteration. 

The  forward  solution  is  re-calculated  using  the  new  updated  values  in  order  to 
determine  the  reconstruction  error, 
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where  nd  is  the  total  number  of  detectors  (56),  6*al  and  /5ccai  are  the  predicted 
detector  data  and  6*hs  and  are  the  simulated  experimental  detector  data  for 

the  excitation  light.  Again,  erg  and  ajAC  are  the  typical  standard  deviations  of  noise  in 
the  measurements  taken  to  be  0.1°  and  1%,  respectively.  For  comparison,  our  multi¬ 
pixel  measurements  have  errors  on  the  order  of  0.2°  and  1.8%.  The  values  of  ag  and 
ajAC  are  scale  factors  that  describe  the  confidence  in  the  measurement.  Every  node 
on  the  grid  will  yield  an  equation  and  the  only  known  variables  are  at  the  detector 
positions.  Since  there  are  many  more  nodes  than  detectors,  the  inversion  scheme  is 
underconstrained  and  not  guaranteed  to  converge  on  the  actual  optical  property  map. 

The  entire  procedure  of  iteratively  adjusting  fxaxf  to  minimize  the  xl  error  contin¬ 
ues  until  a  predetermined  convergence  criterion  is  met.  The  criterion  is  met  if  any  of 
the  following  quantities:  (i)  the  value  of  xl>  (H)  the  absolute  change  in  or  (iii)  the 
relative  change  in  xl  become  lower  than  1.0  xlO-5  or  a  maximum  of  50  iterations  is 
reached.  Typical  times  to  complete  the  inverse  solution  are  approximately  45  minutes 
on  a  Ultra  2  Sun  workstation. 

If  only  a  map  of  /u0i/  is  desired,  the  program  stops.  Otherwise,  the  map  of  fiaxf  is 
used  in  order  to  compute  a  second  map  of  either  <f>  or  r.  Currently,  the  computation 
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of  the  maps  of  j)  and  r  are  conducted  in  two  separate  programs  in  order  to  evaluate 
their  individual  performances. 

After  the  loop  in  which  the  map  / iaxf  is  calculated,  the  forward  solution  is  now 
calculated  to  obtain  6m  and  I™c  measurements  at  the  detector  positions  for  the  emis¬ 
sion  light  using  the  homogeneous  map  described  above  except  with  the  updated  val¬ 
ues  of  //0i/.  Next,  two  Jacobian  matrices  are  constructed:  3 (6m,<p)  and  3(I™C,4>) 
for  the  reconstruction  of  <f>,  or  J(0m,r)  and  3(I%c,t)  for  reconstruction  of  r.  Each 
element  of  these  Jacobian  matrices  is  defined  as  jij  =  36™/3(f>j,  jij  =  dl™c./3<f)j, 
jij  =  36™  j  dr j,  and  jij  =  dl^Q./dTj,  respectively.  Again,  each  element  of  the  Jaco¬ 
bian  matrix  describes  the  response  of  the  source-detector  pair  at  position  i  to  changes 
in  <f>  or  r  at  each  pixel  j.  The  partial  derivatives  are  approximated  as  described 
above  where  36™ /34>  ~  A6m/A(f>,  31™c/d<]>  ~  AI™C/ Aj>,  36m  jdr  ~  A6™  /  At ,  and 
dlZc/dr  «  AIZc/At. 

Equations  3.5  and  3.6  provide  updates  of  (j  and  r  by  adding  A4>  and  At,  respec¬ 
tively,  to  the  values  of  <f>  and  t  from  the  previous  iteration, 
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Again,  the  forward  solution  is  re-calculated  using  the  updated  values  in  order  to 
determine  the  reconstruction  error, 
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where  6™h  I™ceal  are  the  predicted  detector  data  and  6^ls  and  /”co6s  are  the  simulated 
experimental  detector  data  for  the  emission  light.  The  values  for  ag  and  ajAC  are  the 
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same  as  in  the  absorption  reconstruction,  0.1°  and  1.0%,  respectively.  The  entire 
procedure  of  iteratively  adjusting  <j>  or  r  to  minimize  the  xh  error  continues  until  the 
same  convergence  criterion  as  above  is  met  or  50  iterations  have  passed. 
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4.  PART  IV-RECONSTRUCTION  OF  ABSORPTION 
AND  FLUORESCENCE  QUANTUM  EFFICIENCY  AND 
LIFETIME  IN  A  SCATTERING  MEDIUM 

4.1  Results  for  the  Reconstruction  Algorithm 

Three  different  types  of  simulations  were  separately  performed.  They  included 
reconstructions  on  the  basis  of  either  absorption,  Hax},  quantum  efficiency,  or  life¬ 
time,  r.  In  all  three  classes  of  simulated  experiments,  an  optical  heterogeneity  was 
0.063  cm2  which  constituted  0.39%  of  the  total  area.  During  the  reconstructions  of 
both  4>  and  r  unphysically  high  optical  property  values  around  the  periphery  and 
especially  near  the  source  locations  were  obtained.  These  values  were  replaced  with 
the  average  background  value  with  a  constraint  statement  inside  the  program. 

For  all  the  reconstructions,  values  of  jj,axf  range  from  0.010  to  0.200  cm-1.  These 
values  of  absorption  correspond  to  an  ICG  concentration  of  approximately  0.076  to 
1.53  /iM  [4].  These  ICG  concentrations  are  well  below  lethal  levels  and  are  approxi¬ 
mately  5  to  60  times  lower  than  the  therapeutic  concentrations  currently  administered 
with  many  photodynamic  agents  [42,  43]. 

4.2  Reconstruction  of  the  Absorption  Coefficient  from  Excitation  Light 

This  section  describes  the  simulations  performed  on  the  basis  absorption  due  to 
fluorophores,  fiUxf,  at  the  excitation  wavelength.  Figure  4.1  illustrates  both  simulated 
and  reconstructed  data  as  an  optical  tissue  heterogeneity  is  moved  diagonally  towards 
the  detectors  at  three  different  object  locations.  The  heterogeneity  contains  twenty 
fold  more  absorbing  dye  than  its  surroundings  (see  Table  4.1).  The  experiment  shows 
that  the  heterogeneity  location  is  successfully  reconstructed.  However,  since  the  pro¬ 
gram  preserves  the  average  pixel  value,  a  smoothing  of  the  reconstructed  image  is 


34 


observed  causing  the  magnitude  of  ^ai/  for  both  the  object  and  the  surroundings 
to  be  incorrect  as  shown  in  Table  4.1.  As  the  object  moves  closer  to  the  detectors, 
less  error  in  the  reconstruction  is  obtained  since  the  value  of  naxf  inside  the  object 
increases  while  the  value  in  the  background  decreases.  In  fact,  the  value  of  /j.axf  at 
the  closest  object  position  is  0.068  cm-1  (case  c)  is  twice  the  values  of  the  farthest 
position  (0.034  cm-1,  case  a).  This  reduction  of  error  is  expected  since  there  is  more 
signal  yielding  a  larger  disturbance  to  be  measured  at  the  detectors.  In  other  words, 
a  heterogeneity  present  close  to  the  source  and  detector  will  have  more  influence  on 
the  inversion  problem  since  the  pixels  closer  to  the  source  and  detectors  have  more 
influence  on  the  measured  signal  making  the  solution  better  defined.  Successful  re¬ 
constructions  have  also  been  performed  for  a  more  realistic  uptake  of  10:1  as  well  as 
higher  resolutions  (33  x  33  grid)  [44]. 

4.3  Reconstruction  of  Quantum  Efficiency  from  Emission  Light 

The  performance  of  the  inverse  algorithm  for  reconstructing  maps  of  quantum 
efficiency,  <f>,  using  detector  data  at  the  emission  wavelength  is  discussed  in  this  sec¬ 
tion.  For  all  the  inversions  discussed  here,  the  exact  value  of  f, iaxf  inside  the  object 
was  input  into  the  initial  homogeneous  guess.  This  procedure  causes  the  inversion 
algorithm  to  converge  after  one  iteration  in  the  excitation  loop  and  the  exact  value  of 
Haxj  to  be  used  in  the  emission  loop  in  order  to  only  investigate  the  reconstruction  on 
4>.  Table  4.2  lists  all  the  optical  properties  and  experimental  parameters  used  in  the 
forward  solutions  to  generate  the  simulated  detector  data  and  compares  the  values 
from  the  forward  solution  to  those  obtained  by  the  inversion  solution. 

Figure  4.2  shows  the  (a)  actual  solution,  (b)  the  reconstruction,  and  (c)  the  value 
of  <£  inside  the  object  as  a  function  of  iteration  number.  The  experiment  ended  after 
50  iterations  without  meeting  any  of  the  convergence  criteria.  The  optical  property 
map  of  the  reconstruction  shows  that  the  location  of  the  object  is  successfully  located 
but  the  value  of  cf)  inside  the  object  only  reaches  0.082  (18%  smaller  than  the  actual 
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value  of  0.100).  Also,  the  background  pixel  values  are  noisy  with  an  average  value  of 
0.012  which  is  20%  larger  than  the  actual  value  of  0.010. 

4.4  Reconstruction  of  Lifetime  from  Emission  Light 

The  performance  of  the  inversion  algorithm  for  reconstructing  lifetime,  r,  using 
detector  data  at  the  emission  wavelength  is  discussed  in  this  section.  Similar  to  the 
reconstructions  of  <f> ,  all  the  inversions  in  this  section  use  the  exact  value  of  fx axf 
inside  both  the  object  and  background  for  the  initial  optical  property  map.  This  pro¬ 
cedure  causes  the  inversion  algorithm  to  converge  after  one  iteration  in  the  excitation 
loop  and  the  exact  value  of  naxf  to  be  used  in  the  emission  loop  to  investigate  the 
reconstruction  of  r  only. 

The  lifetime  reconstruction  for  a  10:1  uptake  of  dye,  the  current  uptake  ratio  for 
available  contrast  agents,  was  investigated.  Table  4.3  lists  all  the  optical  properties 
and  experimental  parameters  used  in  the  forward  solution  to  generate  the  detector 
data  and  compares  the  values  from  the  forward  solution  to  those  obtained  by  the 
inverse  solution.  This  particular  simulation  had  an  arbitrary  lifetime  value  of  1.00 
ns  inside  the  object  and  a  10.0  ns  everywhere  else.  Figure  4.3  (a)  shows  the  actual 
spatial  map,  (b)  the  reconstructed  spatial  map,  and  (c)  the  average  value  of  r  as  a 
function  of  iteration  number  inside  the  object.  Because  r  inside  the  object  was  lower 
than  in  the  background,  the  object  location  is  defined  as  the  pixel  with  the  lowest 
value  of  r.  Both  the  location  and  magnitude  of  r  inside  the  object  are  correctly  found 
since  the  reconstructed  value  for  pixel  (11,7)  is  0.996  ns,  only  0.04%  lower  than  the 
actual  value.  The  reconstruction  for  the  background  gives  a  good  convergence  of  10.1 
ns  even  though  the  background  contains  a  small  amount  of  noise  which  is  symmetric 
about  the  object.  This  noise  is  most  likely  an  artifact  from  the  large  fluence  that 
emerges  from  the  four  source  locations. 
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4.5  Discussion 

Numerical  simulations  were  conducted  in  order  to  examine  the  resolution  and 
accuracy  of  the  reconstruction  method.  These  results  show  the  potential  of  photon 
migration  for  the  optical  detection  of  a  heterogeneity  inside  a  scattering  medium.  The 
simulated  experiments  based  on  /j,axf  can  locate  an  absorbing  heterogeneity  inside  a 
tissue  simulating  phantom,  however  the  magnitude  of  /j,axf  inside  the  heterogeneity 
is  smaller  than  the  actual  solution.  The  simulated  experiments  also  show  the  ability 
to  reconstruct  the  interior  optical  property  maps  for  both  quantum  efficiency  and 
lifetime  using  detector  information  at  the  emission  wavelength.  Reconstructing  maps 
of  lifetime  not  only  provides  object  location  but  also  may  provide  information  regard¬ 
ing  the  environment  (see  Section  1.16  for  discussion  on  the  Stern- Volmer  equation). 
This  information  about  the  local  bio-chemistry  can  help  to  differentiate  normal  from 
diseased  tissue. 

To  couple  these  inversion  algorithms  with  experimental  measurements,  all  the 
optical  coefficients  of  the  background  need  to  be  scaled  by  a  factor  of  (3/2 )1/2  to 
account  for  the  two  dimensional  nature  of  the  algorithm  and  the  three  dimensional 
nature  of  the  experiment  [45].  The  noise  for  the  fluorescent  measurements  needs  to  be 
determined  since  it  is  unlikely  to  be  the  same  as  the  excitation  data  (0.1°  for  6X  and 
1%  for  I%c)-  Also,  both  inversions  assume  that  the  noise  is  uniform  over  the  entire 
region  which  may  not  be  the  actual  situation,  and  updates  to  the  inversion  scheme, 
which  take  the  noise  as  a  function  of  position  into  account,  need  to  be  implemented. 

To  further  improve  the  reconstructions,  detector  information  at  multiple  mod¬ 
ulation  frequency  needs  to  be  incorporated  into  the  algorithm.  Multiple  frequency 
information  will  provide  more  detector  information  helping  with  the  underconstrained 
nature  of  the  inversion  problem. 
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4.6  Conclusion 

Absorption  and  fluorescent  methods  of  inducing  contrast  were  examined  using 
single  pixel  measurements.  The  results  in  this  study  show  that  fluorescence  has 
improved  contrast  over  absorption  due  to  the  additional  mechanism  from  the  kinetics 
of  the  fluorescence  decay  process.  Also,  it  was  demonstrated  that  lifetime  differences 
alter  the  propagation  of  the  detected  signal.  Since  targeted  delivery  of  a  contrast 
agent  may  not  always  be  feasible,  lifetime  sensitive  dyes  which  are  specific  to  different 
environmental  conditions  found  in  normal  and  diseased  tissue  may  enhance  detection. 
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Figure  1.1.  Dependence  of  absorption  on  wavelength  illustrating  the  “therapeutic 
window”  in  which  the  propagation  of  light  through  tissue  is  high.  Reproduced  from 
Reference  [7]. 


Figure  1.2.  In  the  time  domain,  an  impulse  of  light  is  launched  into  a  scat¬ 
tering  medium.  The  detected  pulse  is  broadened  due  to  the  increase  in  photon 
“times-of-flight .  ” 
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Figure  1.3.  In  the  frequency  domain,  a  sinusoidally  modulated  source  is  launched 
into  a  scattering  medium.  The  detected  signal’s  phase  is  shifted  and  amplitude  is 
decreased  relative  to  the  incident  light. 
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(a)  Detector  ^  Detector 
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Figure  1.4.  Light  propagation  in  a  (a)  non-scattering  and  (b)  scattering  medium. 


Figure  1.5.  Schematic  describing  the  propagation  of  a  photon  density  wave  across  a 
medium  in  the  absence  and  presence  of  a  heterogeneity  with  optical  contrast. 
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Figure  1.6.  Schematic  illustrating  frequency  domain  photon  migration  measurements 
of  a  sinusoidally  modulated  excitation  light  source  at  position  rs  and  the  detected 
excitation  and  emission  (fluorescent)  signals  at  position  when  contrast  is  owing  to 
(a)  absorption,  (b)  perfect  uptake  of  a  fluorescent  dye,  and  (c)  imperfect  uptake  of 
fluorescent  dye. 
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Figure  1.7.  Jablonski  diagram  illustrating  the  electronic  transitions  associated  with 
absorption,  fluorescence  and  phosphorescence. 


COO'  COO 


Figure  1.8.  Molecular  structure  of  heme  showing  the  characteristic  heterocyclic  rings 
that  are  common  to  porphyrin  compounds. 
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Figure  1.9.  Molecular  structure  of  Indocyanine  Green  (ICG). 


Figure  2.1.  Schematic  of  the  cylindrical  tank  used  in  the  single  pixel  measurements  of 
modulation  phase  and  modulation  ratio  for  absorbing  and  fluorescing  heterogeneities. 


Splitter 

Figure  2.2.  Frequency  domain  system  for  single  pixel  measurements  using  a 
Ti:Sapphire  laser  as  the  light  source. 
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Figure  2.3.  Phase  contrast,  A0(0preSence  —  0 absence ),  as  a  function  of  object  position 
for  both  a  perfect  absorbing  object  (open  symbols)  and  a  fluorescent  object  (closed 
symbols)  measured  at  80,  160  and  240  MHz. 


Figure  2.4.  Experimental  measurements  of  (a)  phase  contrast,  A 6  (6 presence— 6 absence), 
and  (b)  modulation  contrast,  AM  (Mpresence / Mat,sence) ,  versus  object  position  for  a 
micromolar  solution  of  DTTCI  at  modulation  frequencies  of  80,  160  and  240  MHz. 
A  9.5  mm  diameter  heterogeneity  was  used  to  hold  the  dye  solution  inside  the  20  cm 
diameter  tank. 
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(a)  (b) 


Figure  2.5.  Experimental  measurements  of  (a)  phase  contrast,  A 9  (0presence  —  9 absence), 
and  (b)  modulation  contrast,  AM  ( Mpresence  /Mab$ence ) ,  versus  object  position  for  a 
micromolar  solution  of  ICG  at  modulation  frequencies  of  80,  160  and  240  MHz.  A 
9.5  mm  diameter  heterogeneity  was  used  to  hold  the  dye  solution  inside  the  20  cm 
diameter  tank. 
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Detector 


Figure  3.1.  Schematic  of  the  two-dimensional  simulated  tissue  phantom  showing  the 
location  of  the  4  sources  and  56  detectors  around  the  perimeter  for  the  17  x  17  mesh 
on  a  16  cm2  area.  The  small  square  object  inside  the  phantom  represents  a  simulated 
diseased  tissue  volume. 
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Figure  3.2.  Flow  diagram  of  the  reconstruction  algorithm.  The  excitation  loop  re¬ 
constructs  Hax)  and  the  emission  loop  either  reconstructs  </>  or  r. 


47 


2  4  ••  10  12  14  16  2  4  6  •  10  12  14  16 

gridnuntwr  gridnuntMr 
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Figure  4.1.  Reconstructed  spatial  maps  of  an  absorbing  heterogeneity  as  a  function 
of  position  for  a  20:1  uptake  ratio  of  dye  inside  the  object.  Figures  (a),  (b)  and  (c)  are 
actual  spatial  maps  and  Figures  (d),  (e)  and  (f)  are  their  corresponding  reconstruc¬ 
tions.  The  actual  values  of  iiax,  in  the  object  and  the  surroundings  are  0.200  cm-1 
and  0.010  cm-1  respectively. 


Figure  4.2.  Reconstructed  spatial  map  of  fluorescent  quantum  efficiency,  4>,  for  a 
10:1  uptake  of  dye  inside  the  heterogeneity.  Figure  (a)  represents  the  actual  image, 
(b)  shows  the  corresponding  reconstruction  and  (c)  shows  the  value  of  <j>  inside  the 
object  as  a  function  of  iteration  number.  The  actual  values  of  <f>  in  the  object  and  the 
surroundings  were  0.100  cm-1  and  0.010  cm-1,  respectively,  and  the  reconstructions 
yielded  values  of  0.082  cm-1  and  0.012  cm-1,  respectively. 
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Figure  4.3.  Reconstructed  spatial  map  of  fluorescent  lifetime,  r,  for  a  10:1  uptake 
of  dye  inside  the  heterogeneity.  Figure  (a)  shows  the  actual  image,  (b)  is  the  cor¬ 
responding  reconstruction  and  (c)  depicts  the  convergence  of  r  in  the  object  as  a 
function  of  iteration  number.  The  actual  lifetimes  for  the  object  and  background  are 
1  and  10  ns  respectively.  The  reconstruction  converges  correctly  to  1  ns  in  the  object. 
The  background  converges  to  10.067  ns  which  is  0.67%  larger  than  the  actual  value 
of  10  ns. 
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Table  4.1  Optical  properties  and  experimental  parameters  used  as  inputs  for  the 
forward  problem  along  with  values  of  fiaxf  obtained  from  the  reconstructions  for  three 
different  object  locations. 


actual 

reconstructed 

background 

object 

background 

ob 

ect 

f^Sx 

(cm-1) 

(cm-1) 

position 

/*«,/ 

(cm  !) 

(cm-1) 

position 

a 

0.010 

10.0 

0.200 

(11,7) 

0.011 

0.034 

(11,7) 

b 

0.010 

10.0 

0.200 

(13,5) 

0.011 

0.042 

(13,5) 

c 

10.0 

0.200 

(15,3) 

0.010 

0.068 

(15,3) 

Table  4.2  Optical  properties  and  experimental  parameters  used  as  inputs  for  the 
forward  problem  along  with  values  of  <j>  obtained  from  the  reconstructions. 


actual 

reconstructed 

background 

object 

background 

object 

Maxi,  M ami 

(cm-1) 

A*«,/ 
(cm  x) 

ftsx  j  f^Sm 

(cm-1) 

(cm-1) 

(ns) 

4> 

Max/ 

(cm-1) 

<t> 

<t> 

4> 

0.000 

0.020 

10.0 

0.002 

1.00 

0.010 

0.200 

0.100 

0.012 

0.082 

Table  4.3  Optical  properties  and  experimental  parameters  used  as  inputs  for  the 
forward  problem  along  with  values  of  lifetime  obtained  from  the  reconstructions. 


actual 

reconstructed 

background 

object 

background 

object 

(cm-1) 

Max/ 
(cm  J) 

ftsx  j  ftsm 
(cm-1) 

a,  mi 

(cm-1) 

d> 

(ns) 

Max/ 

(cm-1) 

(ns) 

(ns) 

(ns) 

0.000 

0.020 

10.0 

0.002 

0.034 

10.0 

0.200 

0.100 

10.1 

0.996 
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Role  of  higher-order  scattering 
in  solutions  to  the  forward  and  inverse 
optical-imaging  problems  in  random  media 

E.  M.  Sevick-Muraca,  D.  L.  Heintzelman,  J.  Lee,  T.  L.  Troy,  and  D.  Y.  Paithankar 


From  analytical  and  numerical  solutions  that  predict  the  scattering  of  diffuse  photon  density  waves  and 
from  experimental  measurements  of  changes  in  phase  shift  B  and  ac  amplitude  demodulation  M  caused 
by  the  presence  of  single  and  double  cylindrical  heterogeneities,  we  show  that  second-  and  higher-order 
perturbations  can  affect  the  prediction  of  the  propagation  characteristics  of  diffuse  photon  density  waves. 
Our  experimental  results  for  perfect  absorbers  in  a  lossless  medium  suggest  that  the  performance  of  fast 
inverse-imaging  algorithms  that  use  first-order  Bom  or  Rytov  approximations  might  have  inherent 
limitations  compared  with  inverse  solutions  that  use  iterative  solutions  of  a  linear  perturbation  equation 
or  numerical  solutions  of  the  diffusion  equation.  ©  1997  Optical  Society  of  America 

Key  words:  Photon  migration  imaging,  image  reconstruction,  frequency-domain,  inverse  problem, 
biomedical  optical  imaging. 


1.  Introduction 

With  the  development  of  near-infrared  (NIR)- 
emitting  laser  diodes  and  the  realization  that  NIR 
light  can  travel  several  centimeters  through  tissue, 
numerous  groups  have  embarked  upon  development 
of  optical  imaging  as  a  new  medical  imaging  modal¬ 
ity.  Although  approaches  such  as  monitoring  the 
vanishingly  small  component  of  coherent  light  with 
sophisticated  techniques  of  Kerr  filtering,1  time  gat¬ 
ing,2  and  conservation  of  light  polarization3  vary, 
other  approaches  focus  on  monitoring  the  predomi¬ 
nant  optical  signal  re-emitted  from  tissues:  the 
multiply  scattered  signal.  Continuous  wave,4  time 
domain,5  and  frequency  domain6"8  measurements  of 
multiply  scattered  light  have  been  performed  in  sim¬ 
ulated  and  experimental  tissue  phantom  studies  as 
well  as  in  human  studies.9  However,  the  use  of 
these  measurements  for  reconstructing  maps  of  in¬ 
ternal  tissue  optical  properties  for  diagnostic  imaging 
has  been  problematic  because  the  geometric  correla¬ 
tion  between  incident  and  detected  radiation  is  de¬ 
stroyed  by  multiple  scattering. 

While  some  investigators  have  used  direct  image 
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reconstruction  approaches  that  use  measured  optical 
data  to  directly  form  an  image, 9-11  others  have 
sought  full  solution  to  the  inverse-imaging  problem, 
which  relates  external  time-dependent  measure¬ 
ments  made  at  the  periphery  to  an  optical  property 
map  of  interior  volumes  through  the  optical  diffusion 
equation.6-912-16  Two  approaches  for  solving  the  in¬ 
verse  solution  have  been  adopted.  In  the  first  ap¬ 
proach,  the  fluence  associated  with  a  propagating 
photon  density  wave  launched  at  source  position  p, 
and  detected  at  the  tissue  periphery  at  position  pd 
[denoted  4>(p„  prf)]  is  related  to  the  fluence  assumed 
in  the  absence  of  any  optical  heterogeneities  [denoted 
‘fyncfp*.  Pd)]  and  the  internal  optical  property  pertur¬ 
bation  map  A(p)  through  a  Fredholm  integral  equa¬ 
tion  of  the  first  kind: 

‘tKp,,  Pd)  =  $inc(pt,  Pd)  +  J  G( p,  Pd)4»o(p«  p)A(p)d3p.  (!) 

where  G(p,  pd)  is  the  Green’s  function  to  the  diffusion 
equation, 

Gip,  pD)  =  4^p-  pD)  exp{t[(-M,  +  iu/cn)/D]m{ p  -  pD)l. 

describing  the  propagation  of  light  from  position  p  to 
the  detector  at  pd,  and  <t>0  is  the  incident  wave  im¬ 
pinging  at  position  p.  If  one  assumes  that  the  inci- 
dent  wave  impinging  upon  the  position  p  can  be- 
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approximated  by  the  fluence  predicted  in  the  absence 
of  heterogeneities  (i.e.,  4>0  -*  ®°m 

approximation),  upon  measurement  of  fluence  at  a 
variety  of  source-detector  separations  and  upon  dis¬ 
cretisation  of  Eq.  (1),  an  optical  map  of  perturbations 
A(p)  can  be  obtained  to  a  first-order  approximation. 
However,  the  Bom  approximation  4\>  ♦ine  does  not 
account  for  second-  and  higher-order  effects  that 
might  arise  from  the  rescattering  of  photon  density 
waves  associated  with  neighboring  inhomogeneities. 
In  addition,  this  approach  assumes  that  perturba¬ 
tions  in  optical  properties  at  a  position  pj  do  not 
influence  the  propagation  of  light  from  position  p^+1 
to  detsdor  pd  as  described  by  the  Green’s  function 

^Investigators  use  the  Bom  iterative  method  (BIM) 
to  account  for  strong  perturbations  and  for  second-  and 
higher-order  scattering  effects  by  using  the  first-order 
approximation  of  optical  property  perturbations  A(p)  to 
iteratively  calculate  <t>0  or  the  fluence  incident  at  po¬ 
sition  p.  Yao  et  al.™  have  shown  that  the  BIM  tends 
to  compensate  for  the  underprediction  of  the  scat¬ 
tering  and  absorption  properties  of  a  single  heter¬ 
ogeneity  in  an  otherwise  uniform  medium  that 
would  occur  when  the  Born  approximation  (or  single 
iteration)  is  used.  In  addition  to  the  BIM,  the  dis¬ 
torted  Bom  iterative  method  (DBIM)  represents  a 
refinement  in  that  it  also  recompiles  the  Green’s 
functions  G(p,  pd)  to  reflect  changes  in  light  propaga¬ 
tion  arid  it  should  speed  convergence.  Yet  to  date 
there  has  been  no  investigation  that  addresses  how 
these  iterative  reconstruction  algorithms  affect  the 
resolution  or,  more  precisely,  how  these  linear  per¬ 
turbation  approaches  can  be  used  to  accurately  image 
closely  positioned,  multiple  heterogeneities  between 
which  nonlinear  second-  and  higher-order  scattering 
effects  exist. 

In  contrast  with  this  inversion  approach.  Jiang  et 
all  a  and  Arridge  et  al.u  have  used  full  numerical  so¬ 
lution  to  the  diffusion  equation  that  describes  the  in¬ 
terdependence  of  voxel  optical  properties  and  their 
contribution  to  the  re-emitted  optical  signal  detected 
from  cw  and  frequency  domain  measurements.  With 
this  approach,  the  second*  and  higher-order  scattering 
arising  from  multiple  inhomogeneities,  which  are  not 
accounted  for  in  the  first  iterations  of  Bom  and  Rytov 
approximations,  are  incorporated  in  the  forward  and 
inverse  solutions.  The  inversion  consists  of  the  relat¬ 
ing  of  spatial  changes  of  optical  properties  on  detected 
frequency  domain  measurements  through  numerical 
solution  of  the  diffusion  equation  and  then  solving  an 
update  to  the  map  of  optical  properties  from  the  dif¬ 
ference  between  the  signals  that  are  measured  and 
those  that  are  predicted  by  the  forward  solution  to  the 
diffusion  equation.  Note  that  convergence  on  the  op¬ 
tical  properties  is  achieved  with  these  numerical  ap¬ 
proaches  and  not  with  the  iterative  Bom  and  Rytov 
approaches.  Nonetheless,  while  the  computational 
investment  of  iterative  but  analytically  based  inver¬ 
sions  is  less  than  that  of  full  numerical  solution  of  the 
diffusion  equation,  the  relative  performance  of  these 


two  inverse-solution  approaches  has  yet  to  oe  evaiu- 
ated. 

For  this  reason  we  embarked  on  a  gtody  to  assess 
the  contributions  of  neighboring  heterogeneities  by 
fifing  experimental,  numerical,  and  analytical  com- 
nutations  of  scattered  photon  densitjrwayes  from 
perfect  light-absorbing  cylinders.  Specifically  we  ex¬ 
perimentally  and  computationally  monitor  the  multi- 
pie  scattering  of  photon  density  waves  between  two 
neighboring  perfect  cylindrical  absorbers  embedded  in 
a  tissue-mimicking  scattering  medium  to  assess  the 
higher-order  perturbation  effects  on  a  re-emitted,  de¬ 
tected  photon  density  wave.  In  the  following  we 
briefly  review  the  theory  of  higher-order  perturbation 
effects  as  predicted  from  the  experimental  frequency 
domain  measurements  as  well  as  from  the  Helmholtz 
equations.  In  addition,  we  present  experimental 
measurements  and  numerical  solutions  of  the  optical 
diffusion  equation  that  show  the  errors,  which  can  be 
significant,  introduced  by  neglecting  second-order  ef¬ 
fects.  These  errors  can  define  the  limits  of  resolu¬ 
tion  for  two  perfect  absorbers  in  the  inverse-solution 
approaches  that  do  not  use  BIM,  DBIM,  or  the  full 
numerical  solution  of  the  optical  diffusion  equation. 
We  comment  on  these  effects  when  contrast  is  caused 
by  mechanisms  other  than  a  perfect  absorber. 

2.  Theoretical  Background  for  Perturbations 
Associate**  with  Diffuse  Photon  Density  Waves 
Our  work  to  assess  the  contributions  of  higher-order 
perturbations  has  been  motivated  by  the  work  of 
Boas  et  al .»*  Their  approach  for  image  reconstruc¬ 
tion  from  diffusely  propagating  photon  density  waves 
uses  analytical  expressions  to  describe  the  complex 
fluence  of  a  propagating  photon  density  wave  <!>(p)  in 
the  presence  of  m  heterogeneities  by  superposition  of 
incident  d>lnc(p)  and  scattered  d>ie.t>n(p)  waves,  from 
body  k  and  of  the  order  of  n: 

*  ro 

<b(p>  -  +  22  *«/&»)•  (2) 

The  second-order  effects  caused  by  the  scattering  of  a 
photon  density  wave  between  two  objects  are  illus¬ 
trated  in  Fig.  1,  from  the  work  of  Boas  and  co-workers, 
and  delineate  the  contribution  of  these  effects  to  the 
detected  photon  density  wave  «Wp).  Second-  and 
higher-order-scattering  effects  arise  from  the  rescat¬ 
tering  of  an  incident  wave  between  the  two  bodies. 
Boas  et  al.  assume  that  second-order-scattering  effects 
(denoted  by  the  dotted  line  in  Fig.  1)  are  negligible 
under  most  circumstances.  We  explore  this  assump¬ 
tion,  using  experimental,  numerical,  and  analytical 
predictions  of  second-order  contributions. 

A.  Analysis  of  Experimental  Measurements  of  Photon 
Density  Waves 

To  measure  second-order  contributions  we  experi¬ 
mentally  measured  the  detected  photon  density  wave 
in  the  presence  of  one  light-absorbing  object  alone, 
4>h;  object  two  alone,  4>w;  and  with  both  objects 
present,  as  a  function  of  object  positions  p}  and 
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Fig.  1.  Schematic  illustrating  the  incident  wave  4>|M  originating 
from  the  source  at  position  p.  (dashed  iineeh  the  first-order  scat¬ 
tered  wave  ♦«. u,"*'  s rising  from  the  first  heterogeneity  (solid 
lines):  and  the  second-order  wave arising  from  reecatter 
of  the  ftrst-oidor  wavs  off  the  second  heterogeneity  (dotted  lines). 

p2  and  separation  between  the  two  objects  P1-P2  (see 
Fig.  2).  From  Eq.  (2),  expressions  can  be  written  for 
measurements  of  <P*t  and  <P*2  made  at  detector  posi¬ 
tion  pd  in  the  presence  of  single  objects: 

^((Pd)  *  ‘l*iiw(Prf) +  <*W*i"’,(Pd).  (3) 

m  ^IsslPd)  +  ^wsUt  *(Pd)» 

where  A1  denotes  the  presence  of  the  first  body  alone 
and  k2  denotes  the  presence  of  the  second  body. 
With  simple  algebraic  manipulation  of  Eqs.  (2H4), 


Fig.  2.  Schematic  detailing  tho  geometry  used  in  the  calculation 
of  scattered  waves  from  analytical  aapreseion.  The  centroid  of  the 
cylinder  is  the  origin,  with  *  denoting  the  length,  angle#  denoting 
the  angle  in  the  plane  containing  the  source  and  the  detector,  and 
r  denoting  the  radial  direction. 

predicted  by  measurements  of  an(* 

<*>»« e(Pd): 

fyuaiPrf)  **  +  (^**a(Prf)  “  ^atiPd)* 

Because  we  report  our  results  in  terms  of  phase  shift 
and  amplitude  demodulation.  Eq.  (6)  can  be  written 
as 


j/Affti  sin  0*i  +  Af*t  sin  0*j  —  0inc\ 

» tan’  ^  w  6m  +  Mu  coa  0|J  -  Af(ne  cos  ej  ’ 

MttjJfH)  -  [(Mk i  cos  eM  +  Mm  cos  eM  -  Mint  cos  *i  J*  +  W*t  sin  ®»i  +  w»»  »“  ®»*  ”  8in  9"^^' 


(7) 

f8) 


an  expression  accounting  for  higher-order- scattering 
effects  can  be  written  for  the  photon  density  wave  in 
the  presence  of  both  objects  <t>fti.*a(prf): 

QkUtfPJ  "  ♦*t(Pd)  +  ♦»*(Pd)  “  *l>to*(Pd) 

+  *hMsUiit"at(Pd)  +  ^wstilif"  *(P<)  ■*■•••• 

(5) 


B.  Analytical  Predictions  of  Phase  Shift  and  Amplitude 
Modulation  Owing  to  Two  Perfect  Absorbing  Cylinder* 
The  complex  fluence  describing  the  propagation  of  a 
photon  density  wave  can  also  be  computed  analyti¬ 
cally  from  the  Helmholtz  equations.  1B-,#  Tho  com¬ 
plex  fluence  propagating  from  a  source  point  at  p, 
through  an  infinite  random  medium  and  received  at 
position  p,  in  tho  absence  of  an  optical  heterogeneity, 
is  given  by 


For  this  study,  it  is  assumed  that  third*  and  higher* 
order  scattered  waves  are  considered  to  have  insig¬ 
nificant  contributions  to  the  measured  fluence  when 
compared  with  those  of  first*  and  second-order  scat* 


A. - 1 

terra 


From  frequency  domain  measurements  of  phase 
shift  e  and  ac  amplitude  M,  values  of  the  complex 
fluence  ♦-  Me~u  can  be  obtained  (1)  in  the  presence 
of  the  first  object  alone,  $*i(Prf);  (2)  the  presence  of 

the  second  object  alone,  ♦*j(p<*);  (3)  in  the  presence  of 
both  objects,  d>*ua(Pa);  *“d  <4>  “  ab*fn“ 
inhomogeneities.  <tWprf).  ^.»«fniord»  effects, 
ie.,  (p*)»  are  neghgible,  from  Eq.  (6)  it 

fbllowstnat  measurements  of  4>*ut(Prf)  should  be 


♦ta«(p)  =  Sms(  pj 


4irZHp  -  94) 


(9) 


where  S,0UTtt  describes  the  phase  and  strength  of 
modulation  of  the  source  located  at  portion  p«;  c*  » 
the  speed  of  light  in  the  medium;  end  D  Ie  the  opticei 
diffusion  coefficient  that  is  governed  by  the  sboorp* 
tion  m  and  the  isotropic  scattering  coeflBoentsoi 
the  continuous  or  homogeneous  medium: 


D  m  l/3(p»>  u»')- 


(10) 


V- 


In  this  study,  we  approximate  a  point-modulated 
source  at  the  periphery  of  a  large  cylinder  as  a  point- 
modulated  source  located  at  the  interface  of  a  semi¬ 
infinite  medium,  hi  this  case  the  fluence  is  written 
in  cylindrical  coordinates  and  at  the  surface  (z  =  0) 
given  by17 


4>me(P<f)  ■ 


S«ar»(p.) 


1 

4ttD 


in  1 

[(p.-p„)2  +  (*-2o)2],/2J 


X 

I 


[(P.-Pd)2  +  (2-2o)2]1/2 

(\  1/2 

[(ft  -  PU)1  +  <*  + 

(Ip.  -  Pd)2  +  (!  +  Z,?]'n  ' 


(U) 


where  z0  is  one  isotropic  scattering  length  (l/p't). 

The  first-order  (n  =  1)  scattered  wave  from  the  £1 
infinitely  long  cylinder  in  an  infinite  medium  is  given 
by1518 


^^["■’(Prf)  =  -'I’mp  i)  f  cos(mfl)cos(pz) 


[f  2  .  (~^a+i^/Cn\ 

1/2  } 

r+(  d  i\ 

*J 

IT  2  .  /~P°  + 

1/2  ’ 

lr  (  D  )J 

Pn 

'  Dxl„'(x)Im(y)  -  D^yJJx)Im'(y) 
DxK„'(x)Im(y)  -  Dki'yKJx)lm'(y ) 


(12) 

where  Dk'  is  the  optical  diffusion  coefficient  within 
the  cylinder  k;  Im  and  Km  are  modified  Bessel  func¬ 
tions;  and  parameters  x  and  y  are  given  as 


x  = 


y  = 


/~Pa  +  »<«>/ 

\  D 

/~Po'  +  tfa)/cn\ 

\  Dk'  I 


1/2 


a*; 


and  a*  is  the  radius  of  cylinder  *1.  Radius  pd,  angle 
0,  and  length  z  denote  the  coordinates  of  the  detector 
relative  to  the  center  of  the  infinite  cylinder  *1  (Fig. 
2).  The  incident  wave  upon  the  infinite  cylinder  <t>inc 
is  computed  from  Eq.  (11).  A  similar  expression  can 
be  written  for  4>,cat  *2n"1- 

We  calculated  second-order-scattering  contribu¬ 
tions  by  evaluating  the  scattered  wave  originating 


from  the  second  object  (k2)  as  the  incident  wave  upon 
the  first  object  (*1).  In  other  words, 


♦acaarU)  =  Z  f  cos(mfi)cos(pz) 

m"l  Jo 

m  ’ 
Pd 
2  1 
Pf 

Dxlm'(x)Im(y)-Dki'yIm(x)Im'(y) 


[p2  +  I 

f-Po  + 

D  /J 

—2  .  1 

('-Po  + 

P  +  1 

D  jj 

DxKm'(x)Im(y)  -  Dki'yKm(x)Im'(y) 


ip, 


(13) 


where  the  incident  wave  upon  cylinder  *1  is  now  the 
first-order  scattered  wave  <h?c«t^2n"1  that  is  com¬ 
puted  from  Eq.  (12).  A  similar  expression  can  be 
written  for  <J>KaU2"''2.  Because  the  optical  proper¬ 
ties  and  the  diameter  of  the  two  cylinders  were  iden¬ 
tical  in  this  study,  we  consider  the  case  in  which  Dt' 
—  D2  and  aj  =  a2. 

With  the  approach  described  above,  the  fluence  de¬ 
tected  at  position  pd  can  be  computed  inclusive  of 
first-order-scattering  effects  (i.e.,  <b(prf)  =  <h,nc(prf)  + 
2*-i  <hgeat,*n’’1(Prf)]>  as  well  as  second-order  scatter¬ 
ing  effects  [i.e.,  <b(pd)  =  4>inc(prf)  + 

<J»,eat^"(Prf)]-  Fro®  final  values  of  complex  fluence, 
the  phase  shift  and  amplitude  demodulation  can  be 
predicted  from  the  simple  relationships: 


0(prf)  =  tan'1 


Im[d>(prf)] 

Re[«>(prf)] 


Af(pd)  =  ({Im[<J>(Prf)]}2  +  {Re[<I>(prf)]}2)I/2. 


(14) 

(15) 


3.  Materials  and  Methods 


A.  Experimental  Measurements  of  Phase  Shift  Owing  to 
Two  Cylindrical  Absorbers 

To  experimentally  measure  MklJt2(pd),  Af*i(p<#), 

^*2(Prf)>  ^inc(prf)>  ®*X>2(Pd)»  (**l(Prf)>  ®*2(Prf)»  and 
0inc(pd),  frequency  domain  measurements  were  made 
with  an  apparatus  that  uses  picosecond  pulsed  light 
at  780  nm  with  an  average  power  of  1.3  W.  Details 
of  the  apparatus  are  described  elsewhere.19  The 
phantom  consisted  of  a  plexiglass  cylinder  (with  a 
16.5-cm  diameter  and  20-cm  height)  filled  with  a 
0.5%  Intralipid  solution  (Kabi  Pharmacia,  Inc.,  Clay¬ 
ton,  N.C.).  As  illustrated  in  Fig.  3,  light  was  deliv¬ 
ered  to  a  peripheral  point  on  the  cylinder  with  a 
1000-pm  fiber  (Model  HCP-M1000T-08,  Spectron 
Specialty  Optics  Co.,  Avon,  Conn.)  and  collected 
through  a  second  fiber  located  4  circumferential  cm 
from  the  incident  source.  Bakelite  plastic  rods  (di¬ 
ameter,  3.175  mm)  were  painted  black  to  provide 
perfectly  light-absorbing  inhomogeneities.  Mea¬ 
surements  of  0(prf)  and  M(prf)  at  80  MHz  were  con¬ 
ducted  as  the  rods  were  moved  in  tandem  along  the 
plane  perpendicular  to  a  line  connecting  the  source 
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measurements. 


and  detector.  Positioning  was  achieved  with  a  mo¬ 
tion  controller  (Model  PMC200-P,  Newport  Corp.,  Ir¬ 
vine.  Calif.)  and  a  motorized  actuator  (Newport 
Model  850B).  The  actuator  position  was  accurate  to 
within  0.0050  mm.  Phase  and  ac  modulation  were 
recorded  as  the  distance  between  the  center  of  the 
first  absorbing  cylinder,  and  the  wall  of  the  phantom 
was  varied  from  0  to  5  cm  in  20  increments.  When 
the  two  objects  were  moved  in  tandem,  their  separa- 
tion  distances  (Pl-p2)  between  the  centers  of  the  two 
perfect  absorbers  were  6,  10,  and  20  mm.  Three 
measurements  were  taken  at  each  position.  Mea¬ 
surements  of  phase  shift  and  ac  amplitude  demodu¬ 
lation  are  reported  relative  to  the  absence  case  or  0inc 
or  Mine. 

B.  Analytical  Prediction  of  Phase  Shift  and  Amplitude 
Modulation  Owing  to  Two  Absorbing  Cylinders 
In  addition  to  the  measurements  of  second-order  in¬ 
teractions  with  the  experimental  approach  described 
above,  predictions  of  interactions  between  absorbers 
were  computed  analytically  with  Eqs.  (11M15)  with 
a  modified  version  of  an  algorithm  written  by  Boas 
that  is  available  through  the  internet.20  Complex 
fluence  incorporating  first-  and  second-order¬ 
scattering  contributions  was  used  to  calculate  phase 
shift  owing  to  two  infinite  cylinders  that  effectively 
act  as  perfect  absorbers  (m  =  2  cm-1,  p..'  7  1°  cm  ' 
in  a  turbid,  semi-infinite  medium  mimicking  our 
phantom  (p*  =  0.02  cm-1,  p.,'  =  10  cm  ).  Zero 
fluence  boundary  conditions  were  assumed  and  used 
in  the  algorithm. 


C.  Numerical  Prediction  of  Phase  Shift  and  Amplitude 
Demodulation  Owing  to  Two  Absorbing  Heterogeneities 
In  addition  to  experimental  measurements  and  ana¬ 
lytical  calculations,  two-dimensional  finite-element 
computations  were  conducted  to  predict  the  phase 
shift  and  amplitude  demodulation  resulting  from 
first-  and  second-order  interactions.  The  computa¬ 
tions  were  performed  on  an  Ultra  2  Sun  Sparc  work¬ 
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station  with  a  matlab  partial  differential  equation 
toolbox.  We  solved  the  frequency  domain  diffusion 
equation  for  light  propagation  to  determine  the  flu¬ 
ence  from  an  80-MHz  sinusoidally  modulated  light 
source.  The  simulated  phantom  was  an  infinite  cyl¬ 
inder  16.5  cm  in  diameter.  The  optical  properties  of 
the  medium  were  set  to  mimic  the  experimental  con¬ 
ditions  of  0.02  cm,-1  for  the  absorption  coefficient  p„ 
and  10  cm-1  for  the  isotropic  scattering  coefficient 
„  >.  We  modeled  the  nearly  perfect  absorbers  in 
these  computations  as  infinite  cylindrical  rods  and 
we  approximated  by  increasing  the  absorption  coef¬ 
ficient  to  100  times  that  of  the  surrounding  medium 
(u,  =  2.0  cm-1).  The  scattering  coefficient  for  the 
object  was  set  to  that  of  the  surroundings  (|is'  =  10 
cm-1)  to  mimic  the  analytical  computations.  The 
phantom  was  discretized  into  66,048  triangular  ele¬ 
ments  that  contained  33,281  nodes.  A  partial  cur¬ 
rent  boundary  condition  was  used  to  approximate 
light  reflection  and  transmission  across  the  boundary 
of  the  phantom  as  would  be  expected  in  the  experi¬ 
mental  measurements.  Although  our  meshing  did 
not  permit  representation  of  the  perfect  absorber  as  a 
volume  excluded  for  light  transport,  we  approxi¬ 
mated  the  heterogeneity  with  a  high  absorption  co¬ 
efficient  because  others  have  shown  that  the 
propagation  characteristics  are  comparable.18  It  is 
noteworthy  that,  although  these  computations  do  not 
exactly  mimic  the  experimental  measurements  de¬ 
scribed  below,  they  nonetheless  adequately  describe 
the  contributions  of  second-order  scattering  effects. 
The  forward  solution  was  obtained  for  each  absorbing 
cylinder  alone  and  then  in  combination.  The  results 
were  analyzed  in  a  manner  similar  to  the  analysis  0 
the  experimental  results  with  Eq.  (7)  and  (8). 

4.  Results  and  Discussion 


^  Experimental  Results 

figures  4(a)-4(c)  and  5(a)-5(c)  illustrate  the  p.iase 
hifl  and  the  amplitude  demodulation  changes  as  a 
unction  of  position  of  the  pair  of  perfect  absorbers 
vhose  centers  are  separated  by  6,  10,  and  20  mm. 
fhe  symbols  denote  individual  measurements  01 
>hase  shift  difference  in  the  presence  of  thetwolud- 
len  objects,  and  the  symbols  connected  by  the  lines 
lenote  the  phase  shift  and  the  amplitude  demodula- 
ion  calculated  from  Eq.  (7)  and  (8)  with  measure¬ 
ments  made  in  the  absence  and  in  the  presence  of  a 
Angle  individual  absorber.  The  measured  values  ot 
jhase  shift  and  amplitude  demodulation  denoted  by 
the  open  symbols  are  therefore  reflective  of  higher- 
irder  contributions,  whereas  the  solid  line  reflects 
jnly  first-order  contributions.  Paired  Students  t 
test  shows  that  there  is  a  significant  difference  (P 
0.005)  between  the  set  of  phase  shift  and  ampli 
demodulation  measurements  and  that  predicted  by 
Eqs.  (7)  and  (8)  for  two  perfect  light-absorbing  objects 
separated  by  6,  10,  and  20  mm  and^^boned  at 
various  distances  from  the  source  and  the  dete 
[Figs.  4(a)— 5(e)X  Our  dso  sho« 

magnitude  of  second-order  effects  is  greatest 
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Distance,  p,  from  source-detector  pair  (cm] 

(a) 


Distance,  p,  from  source-detector  pair  |cm) 

(b) 


(c) 

Fig.  4.  Experimental  values  of  0*i>2“#inc  the  phase  difference  (in 
degrees)  relative  to  the  absence  condition  as  a  function  of  distance 
from  the  source- detector  pair.  The  symbols  denote  individual 
measurements  in  the  presence  of  two  cylinders  separated  by  dis¬ 
tances  of  (a)  6,  (b)  10,  and  (c)  20  mm,  whereas  the  line  connects 
predictions  from  Eq.  (7)  and  measurements  of  0*,,  0*2,  and  0lnc.  The 
error  bars  denote  the  propagation  of  measurement  errors  (standard 
deviation)  associated  with  0*lt  0*2,  and  0inc.  The  x  axis  is  reported 
as  the  distance  between  the  wall  and  the  first  cylinder  (&1).  The 
asterisks  denote  significant  difference  (p  <  0.005,  paired  Student’s 
/-test)  between  the  values  experimentally  measured  and  those  ob¬ 
tained  from  Eq.  (7). 


absorbing  heterogeneities  spaced  6  mm  apart  and 
becomes  smaller  at  10  and  20  mm.  The  paired  f-test 
indicates  that  there  are  significant  differences  at  the 
99.5%  confidence  levels  (indicated  by  the  asterisks  in 


Distance,  p,  from  source-detector  pair  (cm) 
(a) 


Distance,  p,  from  souree-Catactor  pair  (cm) 
(b) 


Distance,  p.  from  source-defector  pair  (cm) 

(C) 

Fig.  5.  Experimental  values  of  MhiJlt/Mint  (in  arbitrary  units! 
relative  to  the  absence  condition  as  a  function  of  distance  from  the 
source-detector  pair.  The  symbols  denote  individual  measure¬ 
ments  in  the  presence  of  two  cylinders  separated  by  distances  of  tat 
6  (b)  10,  and  (c)  20  mm,  whereas  the  line  connects  predictions  from 
Eq.  (8)  and  measurements  olMkl,  Mw,  and  Mlnc.  The  error  bars 
denote  the  propagation  of  measurement  errors  (standard  devia¬ 
tion)  associated  with  Mk „  M»2.  and  Afliw.  The  x  axis  is  reported  as 
the  distance  between  the  wall  and  the  first  cylinder  (fcl). 


Fig.  4)  between  individual  experimental  phase-shift 
values  that  reflect  first-  and  higher-order  scattering 
contributions  to  the  detected  signal  and  those  com¬ 
puted  values  that  are  indicative  of  first-order  effects 
only.  From  the  data  for  6-mm  absorber  separation 
illustrated  in  Fig.  4(a),  it  can  be  seen  that  the  actual 
experimental  measurements  of  phase,  shift  change 
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Olstance.  p,  from  source-detector  pair  (cm) 

(a) 


Fig.  6.  Value*  ofe*1>2-0lM  (in  degrees)  calculated  from  the  ana¬ 
lytical  prediction  of  first-order  (solid  line)  and  inclusive  of  second- 
order  (dashed  Hne)-scattering  effects  as  a  function  of  distance  (cm) 
from  the  source-detector  pair  for  two  absorbing  cylinders  sepa¬ 
rated  by  (a)  6,  (b)  10,  and  (c)  20  mm.  The  *  axis  is  reported  as  the 
distance  between  the  wall  and  the  first  cylinder  (41)  and  the  phase 
shift  reported  relative  to  the  absence  case. 


caused  by  two  objects  ere  smeller  then  those  pre¬ 
dicted  by  Eq.  (7),  in  wbich  second-  and  higher-order 
perturbations  are  not  accounted  for.  At  greater  dis¬ 
tances  from  the  source  and  detector,  agreement  be¬ 
tween  experimental  and  predicted  phase  shift  change 
suggests  that  second-order  effects  may  indeed  be  neg¬ 
ligible,  even  for  the  smallest  separation,  but  only  in  a 
region  where  the  objects’  contributions  to  the  de¬ 
tected  signal  is  comparatively  small.  This  is  reason- 
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Fig.  7.  Values  ofOM.M-flincfin  degrees!  calculated  from  analytical 
prediction  of  first  order  ‘solid  curve)  and  inclusive  of  second  order 
i  dashed  curv  e )  for  two  absorbing  cylinders  separated  by  6  mm  and 
interrogated  at  80,  160.  and  240  MHz. 


able  because  second-order  effects  are  expected  to 
increase  with  proximity  to  the  source  and  the  detec¬ 
tor.  From  our  studies,  it  appears  that  separations 
greater  than  20  mm  are  necessary  for  their  accurate 
resolution  through  perturbative  reconstruction  ap¬ 
proaches  when  perfect  absorbers  are  involved. 

B.  Analytical  Computations 

Our  experimental  results  are  also  validated  by  the 
analytical  predictions  that  account  for  first-  and 
second-order  scattering  contributions.  Figures 
6(a)-6(ci  depict  the  qualitative  trends  seen  in  the 
experimental  data  presented  for  phase  difference  in 
Figs.  4(a )— 4(c).  There  appears  to  be  good  agreement 
between  the  trends  observed  experimentally  with 
various  separation  distances  and  analytical  predic¬ 
tions.  Specifically,  higher-order  contributions  sig¬ 
nificantly  perturb  the  detected  phase  shift  values 
when  the  separations  between  the  center  of  two  ab¬ 
sorbing  cylinders  are  6  mm  [Fig.  6(a)]  and  10  mm 
[Fig.  6(b)].  The  contribution  of  higher-order  scatter¬ 
ing  from  cylindrical  absorbers  is  not  significant  when 
the  separation  distance  is  20  mm  [Fig.  6(c)].  How¬ 
ever  comparison  with  experimental  results  shows 
that  there  are  differences  in  the  absolute  magnitude 
and  shape  of  the  phase  shift  change  versus  °bj 
distance.  These  distances  are  most  likely  caused  by 
the  differences  in  the  boundary  conditions,  geometry, 
and  absorber  strength  between  the  analytical  and  the 
experimental  results.  Nonetheless,  the  trends  con¬ 
firm  experimental  results  that  the  presence  o 
second-order  perturbations  reduces  the  phase  shift 
change  owing  to  two  light-absorbing  bodies.  Ne¬ 
glecting  second-order  effects  could  cause  an  underes¬ 
timation  of  absorption  strength  or  size  when 
reconstruction  algorithms  based  on  first-order  per¬ 
turbations  are  used.  . 

In  addition,  we  investigated  the  variation  oi 
second-order  perturbations  with  modulation  fre* 
quency  as  shown  in  Fig.  7.  The  simulated  data  oi 
phase  shift  change  versus  position  of  the  heterogene¬ 
ities  is  depicted  for  two  cylindrical  absorbers  (d  - 


Dittanci,  p,  from  tourco-datactor  pttr  (cm) 


<b) 


Distance.  p,  from  source -detector  pair  (cm) 


(C) 


(C) 

Fip.  8.  Finite  element  computations  of  A*i>2“einc  <»n  depreesf 
relative  to  an  absence  condition  for  considering  only  first-order 
(solid  line)  and  including  second-order  (dashed  line!  perturbations 
as  a  function  of  distance  (cm)  from  the  source- detector  pair  for 
absorbing  cylinders  separated  by  (a)  6,  (b)  10,  and  (c)  20  mm.  The 
.r  axis  is  reported  as  the  distance  between  the  wall  and  the  first 
cylinder  (*1)  and  the  phase  shift  reported  relative  to  the  absence 


Fig.  9.  Finite  element  computations  of  A^W^nc  arbitrary 
units)  referenced  to  an  absence  condition  for  both  the  first-order 
(solid  line)  and  the  second-order  (dashed  line)  perturbations  as  a 
function  of  distance  (cm)  from  the  source— detector  pair  for  absorb¬ 
ing  cylinders  separated  by  (a)  6.  (b)  10.  and  (c)  20 mm.  The x  axis 
is  reported  as  the  distance  between  the  wall  and  the  first  cylinder 
(fcl)  and  the  modulation  ratio  is  reported  relative  to  the  absence 
case. 


case. 


3.125  mm,  =  2  cm-1)  separated  by  6  mm  in  a 
semi-infinite  medium.  Again,  the  phase  shift  is  re¬ 
ported  relative  to  the  absence  case,  and  the  distance 
is  reported  from  the  wall  to  the  center  of  the  first 
cylinder.  Three  frequencies  were  evaluated:  80, 
160,  and  240  MHz.  The  absolute  value  of  phase 
shift  increases  with  frequency,  but  the  contribution  of 


second-order  effect  remains  roughly  the  same  abso¬ 
lute  value.  Consequently  the  error  in  assuming  neg¬ 
ligible  second-order  scattering  effects  becomes 
smaller  at  increasing  frequencies.  This  is  expected 
because  increased  damping  of  a  rescattered  photon 
density  wave  occurs  at  higher  modulation  frequen¬ 
cies.  While  these  results  intimately  depend  upon 
the  choice  of  optical  properties,  they  nonetheless 
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point  out  that  the  error  in  neglecting  second-order 
scattering  contributions  in  analytically  based  recon¬ 
structions  becomes  smaller  with  increasing  modula¬ 
tion  frequencies.  Of  course  the  error  reduction 
occurs  at  the  expense  of  interrogating  a  smaller  vol¬ 
ume  of  tissue  with  a  source  modulated  at  an  in¬ 
creased  frequency.21 

C.  Finite  Element  Computations 

Our  experimental  results  were  validated  not  only  by 
analytical  predictions  but  also  with  finite  element  com¬ 
putations.  Figures  8(a)— 9(c)  show  the  numerical 
computations  that  correspond  to  the  experimental 
phase  shift  data  presented  in  Figs.  4(a)— 5(c).  Again 
there  appears  to  be  good  agreement  between  the 
trends  of  the  experimental  and  the  analytical  results 
and  the  numerical  solutions.  Also,  higher-order  con¬ 
tributions  perturb  the  detected  phase  shift  values 
when  the  separation  distances  between  the  two  rods 
are  6  mm  [Fig.  8(a)]  and  10  mm  [Fig.  8(b)].  The  con¬ 
tribution  of  higher-order  scattering  from  the  cylindri¬ 
cal  absorbers  is  not  significant  when  the  separation 
distance  is  20  mm  [Fig.  8(c)].  Figures  9fa)-9(c)  show 
the  amplitude  demodulation  relative  to  an  absence 
condition.  The  modulation  data  also  show  results 
similar  to  those  obtained  experimentally.  These  com¬ 
putation  results  confirm  that  the  presence  of  second- 
order  perturbations  is  important  for  two  light¬ 
absorbing  objects  that  are  less  than  20  mm  apart. 

5.  Conclusions 

In  summary,  our  experimental  measurements  show 
that  the  contributions  of  higher-order  scattering  of 
propagating  photon  density  waves  may  not  always  be 
insignificant.  While  our  analytical  and  numerical 
computations  do  not  exactly  reproduce  experimental 
conditions  (i.e.,  two-dimensional  finite  element,  semi¬ 
infinite  geometry,  etc.),  they  nonetheless  demonstrate 
that  the  experimental  trends  can  be  attributed  to 
second-order  effects.  Under  conditions  of  high  con¬ 
trast  caused  by  absorption,  analytical  approaches  to 
the  inverse  imaging  algorithm  may  restrict  the  reso¬ 
lution  and  the  sensitivity  of  biomedical  optical  imaging 
performed  in  the  frequency  domain.  An  analogy  can 
also  be  drawn  for  time  domain  and  cw  reconstruction 
approaches  that  do  not  deploy  full  solutions  to  the 
diffusion  equation  to  account  for  the  interdependence 
of  voxel  optical  properties  on  measured  fluences.  Cer¬ 
tainly  our  results  are  based  on  the  worst-case  scenario 
of  perfect  absorbers  and  may  have  less  impact  on  the 
image  reconstructions  involving  imperfect  absorbers. 
Under  conditions  in  which  multiple  heterogeneities 
are  contrasted  with  their  surroundings  on  the  basis  of 
scattering  or  fluorescence,  the  first-order  perturbation 
assumption  may  not  be  as  restrictive.  Indeed,  if  the 
nonlinearity  associated  with  second-order  scattering 
effects  is  small,  noniterative  Bom  and  Rytov  approxi¬ 
mations  are  especially  attractive  because  image  recon¬ 
struction  is  not  computation  intensive.  Nonetheless, 
our  results  suggest  that  inverse  imaging  algorithms 
that  depend  solely  on  first-order  perturbations  caused 
by  local  changes  in  tissue  absorption  properties  may 
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not  provide  a  reconstruction  as  accurate  as  those  algo¬ 
rithms  that  depend  on  the  full  numerical  calculation  of 
the  diffusion  equation  with  specified  boundary  condi¬ 
tions. 

6.  Appendix  A 


1.  Nomenclature 

ak  Radius  of  cylinder  k  (cm) 
cn  Speed  of  light  through  the  medium  (cm -  si 
D  Optical  diffusion  coefficient  of  homogeneous 
medium  (cm) 

Dk  Diffusion  coefficient  inside  the  cylinder  icmi 

In  Modified  Bessel  function 

In.  Derivative  of  modified  Bessel  function  /„ 

Kn  Modified  Bessel  function 

Kn>  Derivative  of  modified  Bessel  function  Kn 
M  ac  demodulation  of  the  incident  light 
(mW/cm2) 

m  Number  of  objects 

p  Integration  variable  in  Helmholtz  equation  de¬ 

scribing  scatter  from  an  infinite  cylinder 
(cm*1) 

S80urCe  Strength  of  modulated  source  at  position  pv, 
represented  as  a  complex  number  of  amplitude 
and  phase  (mW) 

x  Variable  in  Helmholtz  equation  describing 

scatter  from  an  infinite  cylinder  (cm'1) 
y  Variable  in  Helmholtz  equation  describing 

scatter  from  an  infinite  cylinder  (cm"1) 
z  Axial  direction  of  the  cylindrical  heterogene¬ 

ities  or  length  (cm) 


2.  Greek 

<I>  Photon  fluence  (mW/cm2) 

<I>inc  Photon  fluence  of  incident  wave  or  in  the  ab¬ 
sence  of  heterogeneities  (mW/cm2) 

4>#c<lt  Photon  fluence  arising  from  scattered  wave 
(mW/cm2) 

p  Position  vector  (cm) 

p#  Position  of  the  source  (cm) 

pd  Position  of  detector  (cm) 

p*  Center  position  of  object  k  (cm) 

0  Phase  shift  of  light  wave  (deg  or  rad) 

6  Angle  between  pa  and  pd 

pLa  Absorption  coefficient  of  homogeneous  sur¬ 
roundings  (cm”1) 

Isotropic  scattering  of  homogeneous  surround¬ 
ings  (cm”1) 

|xa»  Absorption  coefficient  of  cylinder  (cm  ) 
Isotropic  scattering  of  cylinder  (cm  l) 


3.  Subscripts  and  Superscripts 

Order  of  perturbation  or  nth-order  scattering  ef¬ 
fect 

h  Index  denoting  cylindrical  object  one,  *  1;  two,  *2; 

or  both  objects,  fcl,  k2 

lnc  Index  denoting  absence  measurement  or  condi¬ 
tion 
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